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Project  Abstract  (Original  Proposal) 

Many  experimental  investigations  have  noted  that  complex  shock  /  boundary  layer 
interactions  exhibit  large-scale,  low  frequency  motion  that  is  distinct  from  that  associated  from 
wall-layer  turbulence.  This  motion  results  in  fluctuating  pressure  and  heat-transfer  loads  that  can 
adversely  affect  the  performance  of  supersonic  missile  systems,  among  others.  Accurate 
predictions  of  such  effects  would  appear  to  require  a  time-dependent  modeling  strategy,  such  as 
large-eddy  simulation  (LES)  or  direct  numerical  simulation  (DNS).  The  application  of  these 
techniques  to  flows  at  practical  Reynolds  numbers  requires  immense  computational  resources, 
and  as  yet,  very  few  studies  have  provided  one-to-one  comparisons  with  experimental  data.  A 
promising  alternative  to  traditional  LES/DNS  are  hybrid  large-eddy  /  Reynolds-averaged  Navier- 
Stokes  (LES/RANS)  models,  which  essentially  function  as  a  large-eddy  simulation  away  from 
solid  surfaces  but  revert  to  a  RANS-based  closure  near  solid  surfaces.  Recent  simulations  of  two- 
and  three-dimensional  shock  /  boundary  layer  interactions  using  a  class  of  “zonal”  LES/RANS 
methods  developed  at  North  Carolina  State  University  have  shown  that  these  techniques  can 
predict  time-averaged  flow  field  properties  as  well  as  or  better  than  RANS  models.  The 
techniques  naturally  capture  both  large-scale  and  small-scale  flow  features,  but  an  open  question 
is  whether  the  dynamics  of  the  shock  /  boundary  layer  interactions  is  predicted  correctly.  This 
tightly  focused  research  effort  will  attempt  to  answer  this  question  by  conducting  detailed  hybrid 
LES/RANS  simulations  of  the  Mach  5,  28  degree  compression-comer  flows  studied  by  David 
Dolling’s  group  at  the  University  of  Texas.  This  data  is  the  most  comprehensive  set  available 
for  studying  shock  interaction  dynamics  and  includes  conditionally-averaged  pitot  and  static 
pressure  distributions,  wavelength,  amplitude,  and  frequency  information  associated  with  the 
shock  structure,  and  cross-correlation  analyses.  Simulations  will  be  conducted  for  this 
interaction  and  will  involve  investigation  of  the  effects  of  variations  in  the  turbulence  modeling 
parameters,  the  extent  of  the  computational  domain,  and  the  degree  of  mesh  refinement  on  the 
predictions.  Data  extraction  will  mimic  that  used  in  the  experimental  measurements  (to  the 
extent  possible).  This  will  enable  the  one-to-one,  time-resolved  comparisons  necessary  for  a 
complete  assessment.  Results  obtained  through  a  successful  validation  should  provide  directions 
for  improvement  in  less-expensive,  lower- fidelity  models  (unsteady  RANS  or  very  large  eddy 
simulation  (VLES))  and  will  perhaps  yield  new  insights  into  the  causes  of  large-scale 
unsteadiness  in  flows  driven  by  shock  /  boundary  layer  interactions. 
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1.  Overview 


This  project  has  resulted  in  the  development  of  a  new  hybrid  large-eddy  /  Reynolds- 
averaged  Navier-Stokes  (LES/RANS)  simulation  strategy  that  models  near-wall  turbulence 
through  unsteady  RANS  concepts.  A  flow-dependent  blending  function,  dependent  on  both 
instantaneous  and  ensemble-averaged  turbulence  quantities,  is  used  to  facilitate  the  RANS  to 
LES  transition,  and  in  contrast  to  other  methods  of  this  type  (e.g.  Detached  Eddy  Simulation), 
the  transition  is  not  directly  dependent  on  the  mesh  spacing.  In  contrast  to  our  earlier  ARO- 
sponsored  work  in  this  scope,  the  new  methodology  does  not  require  the  specification  of  a  model 
calibration  constant.  Rather,  the  use  of  both  ensemble-averaged  and  instantaneous  data  enables 
the  model  to  self-adjust  to  changes  in  the  turbulence  length  scales  -  a  feature  not  found  in  other 
strategies.  The  model  has  been  developed  specifically  for  high-speed  flow  applications  and  has 
been  extensively  tested  for  shock  /  boundary  layer  interactions  as  well  as  scramjet  combustor 
flowfields.  The  model  has  also  been  used  to  predict  turbulent  flow  over  airfoils  near  static  stall. 
A  supplement  to  this  project,  awarded  in  2010,  focused  on  the  development  of  data-mining 
strategies  based  on  LES/RANS  results  and  on  the  use  of  ‘pure’  (non-hybridized)  large-eddy 
simulations  as  a  way  of  identifying  weaknesses  in  the  LES/RANS  closure  assumptions, 
specifically  within  the  ‘transition  region’  where  the  model  shifts  from  (unsteady)  RANS  to  LES. 

This  report  is  organized  as  follows.  Section  2  presents  details  of  the  new  model, 
including  the  rationale  behind  its  development  as  well  as  recent  modifications  designed  to 
improve  its  performance  in  certain  situations.  Sections  3-8  describe  several  validation  studies, 
ranging  from  boundary  layer  simulations  to  fully  3D  shock  /  boundary  layer  interactions  to 
turbulent  flow  over  an  airfoil  near  static  stall.  Section  9  discusses  the  results  of  several  wall- 
resolved  large-eddy  simulations  designed  to  provide  data  to  improve  the  LES/RANS  model,  and 
Section  10  discusses  the  use  of  LES/RANS  data  to  evaluate  RANS  models  for  shock  /  boundary 
layer  interactions.  A  summary  of  the  project  is  presented  in  Section  11.  Publications,  students 
supported,  technology  transfer  activities,  and  other  connections  are  presented  as  the  last  part  of 
this  report. 


2.  Model  Development 

To  anchor  the  description  of  the  new  model,  a  brief  development  of  the  original  NCSIJ 
LES/RANS  model  [1,2],  developed  under  prior  ARO  support,  is  presented  first.  The  original 
model  was  applied  to  a  variety  of  flows  [3-5]  with  good  success,  yet  its  need  for  a  calibration 
step  to  pre-select  a  model  constant  reduces  its  generality. 


2.1.  Original  LES/RANS  model  description 


In  contrast  to  detached  eddy  simulation  (DES)  and  related  strategies,  the  original  approach  is 
designed  to  transition  from  RANS  to  LES  deep  within  the  boundary  layer,  at  approximately  the 
location  where  the  boundary  layer  shifts  from  logarithmic  to  wake-like  behavior.  Menter’s  k-oo 
model,  summarized  in  brief  below,  is  used  as  the  basis. 


d(pk)  d(pkuj ) 
at  dxj 


=  jUtS2  -  P* peak  + 


dXj 


(/z  +  oy//,) 


dk 

dxj 


(2.1.1) 
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dt  8X;  SX  ; 


/  \dco 

OX; 


+  2(l  -  Fx  )pcr(o2 


1  dk  da> 
co  dx,  dx; 


(2.1.2) 


where  S  is  defined  as  the  magnitude  of  the  vorticity  vector  for  the  cases  presented  in  this  report. 
Hybridization  of  the  above  equation  is  accomplished  simply  by  re-defining  the  eddy  viscosity: 

»,=f{  r-+(i-r>„sj,  (2.1.3) 


CO 


with  the  subgrid  viscosity  specified  as 

v,sgs=CMS1,2(q2)1/4A3,2,CM  =0.06  (2.1.4) 

An  estimate  of  the  subgrid  kinetic  energy  is  obtained  by  test- filtering  the  resolved-scale  velocity 
data: 


2  l  ~  ^  - 

q  =-(uk-ukY 


(2.1.5) 


The  response  of  the  model  is  dictated  by  the  blending  function  T ,  a  time-dependent  quantity  that 
reaches  a  value  of  unity  where  an  unsteady  RANS  response  is  desired  and  zero  where  an  LES 
response  is  desired.  It  should  be  mentioned  that,  because  of  the  effect  of  unsteady  strain  rates  in 
enhancing  the  production  of  k  and  00,  one  cannot  expect  that  the  ensemble-averaged  values  of  k 
and  co  will  be  equivalent  to  those  obtained  under  steady-state  assumptions  (even  with  T  equal  to 
1).  In  fact,  the  production  of  k  is  greatly  enhanced  by  unsteady  strain  rates;  the  production  of  co 
less  so. 


The  blending  function  used  in  Refs.  [1]  and  [2]  is  based  on  the  ratio  of  the  wall  distance  d  to  a 
modeled  form  of  the  Taylor  micro-scale: 


f  \ 


r  =  y  1  -  tanh[  5(— ^=7/ 2  - 1)  -  ^] 

V  vS*  J 

9 

(2.1.6) 

II 

(2.1.7) 

where  the  Taylor  micro-scale  is  defined  as 

=  ^jV  1^^ 

(2.1.8) 

-1  K  2 

The  constant  6  =  tanh  (0.98)  shifts  the  balancing  position  (where  , _ ri  =1)  from  r=0.5  to 

VC 

T  =  0.99.  It  should  be  mentioned  that  the  shifting  parameter  cj>  actually  corrects  for  the  fact  that 
the  time-averaged  value  for  co  is  larger  than  the  RANS  value  because  of  the  enhancement  of 
turbulence  production  mentioned  above. 

The  constant  cqis  chosen  to  force  the  average  LES  to  RANS  transition  position  ( T  =  0.99)  for 
equilibrium  boundary  layers  to  occur  at  the  point  where  the  wake  law  starts  to  deviate  from  the 
log  law.  To  determine  ax  for  a  particular  inflow  boundary  layer,  the  following  method  is  used. 
First,  a  prediction  of  the  equilibrium  boundary  layer  is  obtained  (given  free-stream  properties,  a 
specified  wall  thermal  condition,  and  a  value  for  the  boundary  layer  thickness)  from  Coles’  Law 
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of  the  Wall  /  Wake  along  with  the  Van  Driest  transformation.  An  initial  estimate  for  the  outer 
extent  of  the  log  layer  is  defined  by  finding  the  value  of  d+  such  that 


Tln(0  +  C|/ 

(u  \ 

uvd 

lUr  J 

=  0.98 


(2.1.9) 


The  value  of  d+  =  u rd  /  v  that  corresponds  to  this  value  of  d+w  is  then  found  through  the  use  of 
Walz’s  formula  for  the  static  temperature  distribution  within  the  boundary  layer: 


T  T^+(^~TJ_u 
T„  T„  T  u„ 


(r- !) 


f  V 


M: 


(2.1.10) 


The  model  constant  is  then  found  by  the  equivalence  d+  =ax  ,  which  arises  from  the  use  of 
inner- layer  scaling  arguments  for  k  and  co.  Specifically,  in  the  logarithmic  region,  one  has 
co  =  nT  /(/cd^jc^)  ,  where  vis  the  Von  Karman  constant,  Cfl  =  ff  =0.09,  and  uT  is  the  friction 
velocity.  Substituting  this  into  Eq.  2.1.8  and  Eq.  2.17  in  succession,  one  finds  that 

K  T  K 


at  the  balancing  position.  Figure  2.1.1  shows  an  example  of  the  positioning  of  the  blending 


CMd  co 


afv 


M  J  -i 
2  2  1 
ax  v  a ! 


(2.1.11) 


Figure  2.1.1:  Blending  function  and  velocity  profile  in  wall  coordinates 

function  within  a  compressible  boundary  layer.  This  procedure  requires  free-stream  information 
and  a  target  boundary  layer  thickness,  both  of  which  are  specific  to  a  particular  boundary  layer. 
As  such,  the  calibration  procedure  is  not  universal  and  the  constant  that  results  must  be  adjusted 
on  a  case-by-case  basis.  Nevertheless,  Ref.  [2]  in  particular  shows  that  the  model  as  described 
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above  accurately  captures  mean  and  second-moment  compressible  boundary-layer  statistics  over 
a  wide  range  of  Reynolds  numbers,  provided  that  the  outer  layer  is  resolved  adequately  (15-20 
cells  /  boundary  layer  thickness  in  all  directions). 


2.2.  New  LES/RANS  model  development 

Our  development  of  a  new  LES/RANS  model  centers  on  the  construction  of  an  alternative  form 
for  the  blending  function  that  does  not  need  a  problem-specific  calibration.  There  are  problems 
with  using  the  Taylor  microscale  as  a  representative  turbulence  length  scale  far  away  from  the 
wall.  At  high  Reynolds  numbers,  outer-layer  statistics  (when  scaled  properly)  should  be 
independent  of  Reynolds  number  and  thus  of  the  molecular  viscosity.  What  is  needed  is  a  local 
estimate  of  the  outer-layer  integral  scale  (proportional  to  the  boundary  layer  thickness). 
Equating  this  value  with  the  inner  layer  integral  scale  (proportional  to  the  wall  distance)  will  give 
a  reasonable  estimate  of  the  location  of  the  outer  part  of  the  logarithmic  region,  and  by  analogy, 
a  location  where  one  might  wish  to  shift  the  model  response  from  RANS  to  LES.  The 
difficulty  is  that  the  LES  model  locally  supplies  only  inner-layer  length  scales  (the  wall  distance 
+  viscous  length  scales)  and  the  grid  scale. 


The  turbulence  transport  equations  may  provide  some  outer-layer  length-scale  information.  As 
one  example,  the  length  scale  that  is  naturally  determined  from  a  k-oi  model  is 


■  1/2 


. k-w 


/ — »1/ 4 

CO 


(2.2.1) 


This  length  scale  approaches  the  inner  scale  Kd  in  the  logarithmic  region  and  decays  toward  zero 
in  the  outer  part  of  the  boundary  layer  and  in  the  viscous  sublayer.  Taking  the  ratio  of  lk_m  to  the 


inner  scale  and  adding  in  a  viscous  velocity  scale  proportional  to  vco  gives  a  potential  form  for 
the  argument  of  a  new  blending  function: 


Vl0v<y  +  ~k 


(2.2.2) 


with  the  constant  (10)  being  chosen  small  enough  so  that  the  behavior  of  X  in  the  logarithmic 
and  outer  layers  is  not  significantly  influenced  by  the  viscosity.  The  length-scale  ratio  X  should 
approach  one  and  larger  in  the  inner  layer  and  should  decay  toward  zero  in  the  outer  part  of  the 
boundary  layer.  The  Menter  k-oo  model  itself  uses  a  ratio  similar  to  X  to  define  its  Fx  function. 
The  delayed  detached-eddy  simulation  (DDES)  model  [6]  also  uses  a  similar  function  to  separate 
attached  boundary  layers  from  free-stream  regions,  and  some  of  our  earlier  attempts  at 
LES/RANS  model  development  also  used  this  basic  form.  An  advantage  of  this  form  is  that 
there  is  a  good  correlation  between  a  particular  value  of  X  and  the  50%  value  of  T  as  determined 
through  the  original  calibration  procedure.  Figure  2.2.1  shows  this  for  the  three  boundary  layers 
considered  in  Ref.  [2],  The  length  scale  ratio  in  this  figure  isLr  =  /l2 .  These  results  show  that  the 


50%  value  of  the  blending  function  correlates  approximately  with  a  value  of  Lr  =  0.6  and  is 
nearly  independent  of  Reynolds  number. 
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There  are  several  problems,  however,  with  the  direct  use  of  A ,  evaluated  using  instantaneous 
values  for  v ,  co,  and  k .  First,  as  shown  in  Figure  2.2.1,  the  decay  rate  of  the  function  in  the 
outer  part  of  the  boundary  layer  is  rather  slow.  If  the  goal  is  to  transition  at  the  outer  edge  of  the 
logarithmic  region  so  that  most  of  the  boundary  layer  is  modeled  as  a  LES,  then  the  functional 


Figure  2.2. 1 :  Correlation  between  blending  function  and  length-scale 
ratio  for  several  flat  plate  boundary  layers 

form  must  be  sharpened  in  some  manner  so  that  the  decay  rate  is  increased.  Forms  for  the 
blending  function  that  may  be  considered  include 


T  =  ^-(l  +  tanh[  C(Ak  - £>)])  or 

(2.2.3) 

r  =  ifi-tanh[C(-L-z>)]' 

z  v  ^  y 

(2.2.4) 

where  C  is  the  sharpening  factor  (normally  >  5) ,  D  is  the  value  for  which  T  =0.5  (the  RANS-to- 
LES  transition  ‘point’) ,  and  m  is  a  power. 

The  other  problems  are  more  subtle.  Because  the  solutions  for  A:  and  co,  as  driven  by  an 
unsteady  velocity  field,  are  intrinsically  noisy,  local  values  for  the  length-scale  ratio  may  differ 
substantially  from  the  time-averaged  values,  which  themselves  may  differ  significantly  from 
those  expected  from  a  RANS-based  analysis  (as  used  in  Figures  2.1.1  and  2.2.1).  There  is  also  a 
strong  coupling  between  the  length  scale  ratio,  the  blending  function,  and  the  transport  equations 
themselves.  The  functional  form  for  the  eddy  viscosity  (Eq.  2.1.3)  depends  on  T,  which  itself 
depends  on  k  and  co .  The  presence  of  T  in  the  turbulence  kinetic  energy  transport  equation,  in 
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particular,  induces  a  large  damping  effect  for  turbulence  production,  which  causes  the  solution 
for  k  to  be  smaller  than  the  RANS  value  (even  where  T  may  be  one  locally).  This  can  be  easily 
seen  by  writing  the  production  term  as 

M,S2  =  [Tp  —  +(\  -  T)pv,  ]V-  =  p-S2  -  (1  -  T)p(--v,  )S2  (2.2.5) 

CO  CO  CO 

A  consequence  of  this  strong  coupling  and  the  diminishing  of  k  near  the  RANS/LES  interface 
is  that  the  fixed-point  solution  of  the  LES/RANS  equations  (determined  by  the  position  of  the 
time-averaged  blending  function  within  the  boundary  layer)  may  depend  on  the  initial 
distribution  of  k .  If,  for  example,  a  RANS  solution  is  used  to  initialize  k ,  the  most  probable 
response  will  be  that  the  blending  function  stays  near  one  throughout  much  of  the  boundary 
layer,  meaning  that  turbulent  fluctuations  will  be  attenuated  significantly.  This  is  the 
operational  mode  required  in  DDES  and  related  techniques,  which  seek  to  treat  the  entirety  of  an 
attached  turbulent  boundary  layer  as  (unsteady)  RANS.  If  the  initial  distribution  for  k  is  scaled 
so  that  the  outer-layer  value  is  near  zero,  then  theT=  0.5  position  blending  function  may  stay 
initially  near  its  target  location,  but  will  likely  move  toward  the  wall  as  the  resolved  turbulence 
intensity  in  the  outer  layer  increases. 

The  original  blending-function  model  described  in  the  previous  section  avoids  some  of  these 
difficulties  in  that  it  is  a  function  of  co  and  not  k ,  the  former  only  being  sensitive  to  the  eddy 
viscosity  through  the  diffusion  terms  and  to  k  through  the  cross-derivative  term.  Thus,  the 
model  yields  stable  fixed-point  solutions  for  the  time-averaged  position  of  the  blending  function. 
However,  because  of  the  absence  of  outer-layer  length  scale  information,  the  original  model 
cannot  be  self-adjusting  and  requires  problem-specific  calibration. 


In  this  work,  we  propose  a  definition  of  the  length  scale  ratio  that  avoids  some  of  the  above 
issues.  Instead  of  only  using  the  modeled  A:  to  define  the  ratio,  we  include  the  resolved 
turbulence  kinetic  energy  kR ,  which  is  obtained  by  ensemble  averaging  the  unsteady  flow  field 
according  to 


pkR-  (pukuk  pi'kpllk ) 

2  p 

The  length  scale  ratio  that  results  is 

(2.2.6) 

_  yj\0vco  +  k  +  kR 

~  ^ N  ^1/4  7  ’ 

Cu  mco 

(2.2.7) 

with  CN  being  a  model  constant.  The  inclusion  of  the  resolved  turbulence  kinetic  energy  ensures 


that  part  of  the  outer-layer  length  scale  estimate  is  more  independent  of  the  positioning  of  the 
blending  function  within  the  boundary  layer.  An  additional  modification  can  be  made  to  reduce 
the  noise  in  the  length-scale  ratio  to  levels  found  in  the  original  model: 


4  —  4  • 


io+k+k* 


4v 


vco  Cl'AKdJco 


(2.2.8) 
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I —  k  +  k 

Here,  we  factor  out  the  quantity  V vao  and  evaluate  the  component  ,  10  h - =r^-  using  ensemble- 

V  vco 

averaged  data  for  co  and  k .  The  new  blending  function  uses  Eq.  2.2.8  in  conjunction  with  Eq. 
2.2.4  (with  C=  15,  D=  1,  and  m  =  2): 

f  ,  A 

(2.2.9) 


1 


1  -  tanh[  15(— j-- 1)] 


In  this  form,  it  is  easy  to  compare  the  new  blending  function  with  the  original.  One  can  re-write 
the  original  form  into  that  of  Eq.  2.2.4: 


1 


To  =-[l-tanh({5  +  tanh  ‘(0.98)}( 


K 


=  — [l-tanh(7.30( 


K 


5  +  tanh_1(0.98)  C'/2(a,)2 
f  d\2 


d_\ 

Xt 


-1)] 


(2.2.10) 


1.46C‘/2(a1) 


\  Xt  . 


■1)] 


Introducing  the  Taylor  microscale  into  Eq.  2.2.10,  one  can  write  the  new  blending  function  as 


fV=i[l-tanh(15( 


1 


K 


f  j  \ 


C/  (10  +  ^±^)c» 


1/2 


V  J 


1)] 


(2.2.11) 


cov 


Equating  the  denominators  in  the  arguments  of  the  original  and  new  blending  function,  it  is  easy 
to  see  that 


ai2._SMio+*i±l*>) 

1.46/r  cov 


(2.2.12) 


In  essence,  the  new  model  replaces  the  problem-specific  constant  ax  with  a  form  that  varies  in 
space  and  also  in  time,  as  the  ensemble  averages  may  not  be  time-invariant  and  as  the  kinematic 
viscosity  is  a  local  quantity.  The  model  constant  CN  must  be  obtained  through  numerical 

experiments.  Considering  the  same  three  compressible  boundary  layers  as  in  [2]  and  comparing 
mean-  and  second-moment  statistics  with  those  obtained  using  the  original  model  (results  shown 
later),  we  have  determined  that  a  value  of  CAr=1.5  gives  acceptable  results.  Figure  2.2.3  plots 

f  „  -  \  1/2 


the  quantity 


C 


_(io 

1  A6kccx  co  v 


N 


versus  the  normalized  wall  distance  for  each  of  the 


^  X  •  I  V/#V  tA/  V  J 

boundary  layers  considered  in  [2],  The  values  of  ax  for  the  Elena  and  Lacharme  [7],  Luker,  et 
al.[8],  and  Smits  and  Muck  [9]  experiments  are  14.16,  23.22,  and  75.24,  respectively.  The  peak 
values  for  each  distribution  are  close  to  unity  but  are  lower  than  unity  where  the  RANS-to-LES 
transition  takes  place  (y  /  8  ~  0.2 ).  Also  shown  in  Figure  2.2.3  are  the  time-averaged  blending 
functions  for  the  new  model  and  for  the  original  model.  In  general,  the  Y  =  0.5  values  are  closer 
to  the  wall  for  the  new  model,  but  the  transition  zone  itself  (0  <  Y  <  1)  is  broader. 
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Figure  2.2.3:  ratio  of  model  ‘constants’  and  blending  function 
distributions  versus  normalized  wall  distance  for  three  compressible 
boundary  layers 


2.3.  Model  modifications 


Since  the  original  development  of  the  new  model,  there  have  been  some  modifications  made  to 
allow  some  control  over  the  RANS-to-LES  transition.  In  some  problems,  there  is  a  need  to 
guarantee  operation  in  RANS-mode  for  various  reasons,  such  as  the  lack  of  sufficient  mesh 
resolution  for  resolving  fluctuations  or  a  user’s  choice  to  not  attempt  to  resolve  fluctuations  in 
certain  regions  of  the  flow.  The  calibration  of  the  current  LES/RANS  model  assumes  that 
sufficient  resolution  is  present  and  that  only  the  inner  20%  or  so  of  a  boundary  layer  will  require 
some  level  of  RANS  modeling.  The  model  therefore  does  not  work  properly  in  situations  where 
a.)  resolved  turbulence  is  not  sustained  through  some  procedure,  or  b.)  the  mesh  resolution  isn’t 
adequate  to  support  resolved  turbulence.  Two  methods  for  controlling  the  transition  have  been 
developed  and  are  described  in  this  section.  The  first  method  modifies  the  blending  function 
argument  as  follows: 

f  11^  d 

1  -  tanh[15(— - -)]  ,  (3  =  max[y90 ,  max*  (Ck  min(l,  — )] 

d,r 


r  =  I 
2 


4  A 


(2.3.1) 


Here,  /?  is  a  function  of  distance  functions  dk  as  calculated  from  a  user-specified  number  ( k  )  of 


solid  surfaces  and  the  global  distance  function  d .  The  user  determines  beforehand  whether  a 
surface  is  to  be  treated  as  a  ‘RANS  wall’  (Ck  =  0)  or  as  a  ‘LES  wall’  ( Ck  =  1).  This  function 


approaches  one  for  cells  near  a  ‘LES  wall’  and  approaches  a  threshold  value  /?0~0.05  for  cells 
near  a  ‘RANS  wall’.  The  effect  is  to  shift  the  RANS-to-LES  transition  location  toward  the  outer 
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edge  of  a  boundary  layer  when  one  does  not  wish  to  resolve  turbulence  within  that  boundary 
layer. 


The  second  method  determines  adequacy  of  the  mesh  resolution  by  comparing  the  estimated 
outer-layer  length  scale  with  the  maximum  mesh  scale  according  to  the  following: 


§  Pouter  Pouter  /  C 


g(l  outer )  =  mint  A » max(l,  D, 


A„ 


l. 


where 


l  =  C 

1  outer  ■ 


1 1 0  vco  +  k  +  k 
C^CQCQ 


=  xd 


£>] 

(2.3.2) 

V  CO 

(2.3.3) 

The  model  constants  Dx  and  D2  are  assigned  values  of  10  and  0.5,  based  on  calibrations  for  flat- 
plate  boundary  layers.  The  mesh  scale 
A  maxis  taken  to  be  the  maximum  spacing 


over  all  three  coordinate  directions.  This 
form  serves  to  shift  the  closure  to  unsteady 
RANS  when  there  is  no  possibility  of 
resolving  the  largest  turbulence  length 
scales.  If  the  maximum  mesh  scale  is 
selected  as  the  outer-layer  scale 
consistently,  then  the  LES/RANS  model 
behaves  similarly  to  detached-eddy 
simulation,  serving  primarily  to  isolate 
turbulent  boundary  layers  from  massively- 
separated  regions.  Figure  2.3.1  shows  the 
effective  RANS  to  LES  transition  position 
as  a  function  of  the  maximum  mesh 
spacing  for  a  flat-plate  boundary  layer.  As 
the  mesh  spacing  gets  larger  than  -10%  of 
the  boundary  layer  thickness,  the  transition 
location  moves  outward  and  the  model 
shifts  more  to  unsteady  RANS. 


Figure  2.3.1:  Length-scale  ratio  versus  Y/8  (Elena 
and  Lacharme  flat  plate  boundary  layer) 


2.4.  Ensemble  averaging 


The  new  LES/RANS  model  requires  ensemble-averaging  of  unsteady  data  to  form  the  turbulent 
length-scale  ratio  required  for  the  blending  function.  These  ensemble  averages  are  currently 
computed  using  an  exponentially-weighted  moving  average:  Qn  =  (\-A)Qn~}  +  AQn 


with  A  =  At/r.  The  time  scale  r  is  defined  as  follows  for  each  of  the  variants  considered  in  this 
study: 


1.  t  =  t 

2.  T  =  mm(t,tres),  t<4tres 


(2.4.1) 


3.  T  =  mm(t,t,.J 
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where  rres  =  L!  ur  is  defined  in  terms  of  the  length  of  the  domain  L  and  the  free-stream  velocity 
ux  Method  3  is  an  exponentially- weighted  moving  average,  with  the  ‘window’  set  as  one  flow¬ 
through  time,  while  Method  1  is  a  conventional  time  average.  Method  2  is  a  blend  of  a 
‘windowed’  average  and  a  conventional  one.  Results  illustrating  the  effect  of  these  variations 
are  described  later  in  this  document. 

3.  Model  Assessment:  Flat  Plate  Boundary  Layers 

3.1.  Baseline  predictions 

Three  supersonic  flat-plate  boundary  layers  experiments  were  simulated  to  calibrate  the  new 
model  and  to  assess  the  sensitivity  of  the  model  to  the  model  constants  and  to  the  method  of 
ensemble-averaging.  Figures  3. 1.1 -3. 1.3  compare  predicted  mean  velocity,  Reynolds  shear 
stress,  and  rms  normal  and  axial  fluctuation  intensities  with  experimental  data  [10-12],  The 
constant  CN  is  set  to  1.5  for  the  new  model,  while  the  constant  Cs  is  chosen  as  15  and  as  5. 

Reynolds  numbers  based  on  the  boundary  layer  thickness  are  5.59e4,  1.78e5,  and  1.58e6  for  the 
Elena  and  Lacharme  [10],  Luker,  et  al.  [11],  and  Smits  and  Muck  [12]  experiments,  respectively. 
The  meshes  each  contain  5.12  million  cells  and  are  designed  so  that  approximately  20  cells  / 
boundary  layer  thickness  are  present  in  the  wall-transverse  directions  (x  and  z).  Significantly 
more  resolution  (>100  cells  /  boundary  layer  thickness)  is  present  in  the  wall  normal  direction. 
Calculations  using  the  original  model  [7]  are  also  presented.  These  results  show  several  trends. 
First,  the  mean  velocity  is  relatively  unaffected  by  the  model  choice.  The  Reynolds  shear  stress 
distributions,  presented  in  terms  of  their  modeled  and  resolved  components,  show  that  the  decay 
in  resolved  Reynolds  shear  stress  is  compensated  for  by  an  increase  in  the  modeled  component. 
The  50%  value  in  the  normalized  modeled  Reynolds  stress  corresponds  approximately  to  the 
50%  value  in  the  time-averaged  blending  function.  With  the  exception  of  the  Smits  and  Muck 
case,  the  differences  among  Reynolds-shear  stress  profdes  predicted  by  each  model  are  small. 
Lowering  the  sharpening  factor  from  15  to  5  leads  to  an  increase  in  the  modeled  Reynolds  shear- 


Figure  3.1.1:  Mean  velocity,  Reynolds  shear  stress,  rms  axial  velocity,  and  rms  normal  velocity 
versus  normalized  wall  distance  (Luker,  et  al.  experiment) 

stress  component,  implying  that  the  position  of  the  time-averaged  RANS-to-LES  transition  shifts 
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further  away  from  the  wall.  All  models  under-predict  the  peak  in  rms  axial  velocity  indicated  in 
all  datasets  at  y/8o  <0.2.  This  appears  to  be  a  consequence  of  the  damping  effect  of  the  RANS 
component  of  the  model.  It  should  be  noted  that  the  modeled  contributions  to  the  Reynolds 
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-  New  model  (GN=1 .5,  Gs=1 5,  Modeled) 

-  New  model  (GN=1 .5,  Gs=5, Resolved) 

-  New  model  (CN=1 .5,  (Cs=  5,  Modeled) 


Figure  3.1.2:  Mean  velocity,  Reynolds  shear  stress,  rms  axial  velocity,  and  rms  normal  velocity 
versus  normalized  wall  distance  (Elena  and  Lacharme  experiment) 


normal  stresses  are  not  included  in  the  rms  velocity  comparisons,  as  the  Boussinesq 
approximation  is  known  to  be  invalid  for  these  components  in  the  near-wall  region.  With  the 
exception  of  the  Smits  and  Muck  case,  the  model  predictions  for  the  rms  normal  velocity  are  in 
good  agreement  with  experimental  data.  In  all  cases,  the  predictions  provided  by  the  new  model 
are  comparable  to  those  of  the  original  model. 


Figure  3.1.3:  Mean  velocity,  Reynolds  shear  stress,  rms  axial  velocity,  and  rms  normal  velocity 
versus  normalized  wall  distance  (Smits  and  Muck  experiment) 


3.2.  Sensitivity  to  CN 

Figure  3.2.1  addresses  the  sensitivity  of  the  predictions  to  variations  in  the  model  constant  CN  . 
The  experiment  of  Luker,  et  al.  is  considered  in  this  evaluation.  As  expected,  a  trend  of 
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increasing  the  modeled  Reynolds  shear-stress  component  and  decreasing  the  resolved  component 
is  evidenced  with  an  increase  in  the  model  constant.  As  the  model  constant  reduces,  more  of  the 
resolved  turbulence  energy  near  the  wall  is  captured,  and  the  prediction  of  rms  axial  velocity,  in 
particular,  improves.  This  effect  is  accompanied,  however,  by  an  increase  in  the  fullness  of  the 
velocity  profile  in  the  outer  part  of  the  logarithmic  region  (0.1  <  y/  5o  <  0.25)  and,  as  shown 
later,  a  decrease  in  the  wall  shear  stress.  The  apparent  over-prediction  of  the  resolved  Reynolds 
shear  stress  for  CN=  1.2  is  partly  a  result  of  this  reduction  in  the  wall  stress.  Even  at  the  highest 

value  of  CN  considered  (2.2),  a  significant  amount  of  fluctuation  energy  is  still  sustained  in  the 

outer  layer.  A  shift  to  a  RANS-like  response  would  require  even  larger  values  of  the  model 
constant. 


Figure  3.2.1:  Effect  of  CN  (Luker,  et  al.  experiment) 


3.3.  Sensitivity  to  ensemble-averaging  method  and  Cs 

The  sensitivity  of  the  predictions  to  the  method  of  ensemble-averaging  and  to  the  selection  of 
Cs  is  illustrated  in  Figures  3.3.1a  and  3.3.1b,  corresponding  to  the  Luker,  et  al.  and  Smits  and 

Muck  experiments.  Comparisons  with  experimental  data  are  again  presented  for  the  resolved 
Reynolds  shear  stress,  but  the  velocity  profiles  are  shown  in  inner-layer  coordinates  using  the 
Van  Driest  I  transformation.  This  transformation  highlights  the  smoothness  of  the  velocity 
profile  in  the  inner  region  and  indicates  the  degree  to  which  a  Tog-law  mismatch’  is  observed. 
The  sensitivity  of  the  predictions  to  the  method  of  ensemble-averaging  is  not  high,  but  the 
conventional  average  (“Weighting  #1”)  places  the  time-averaged  RANS-to-LES  transition 
somewhat  further  away  from  the  wall.  This  implies  that  some  of  effects  of  the  initial  distribution 
of  turbulence  variables  are  maintained  when  conventional  averaging  is  used.  Weightings  #2  and 
#3  provide  very  similar  solutions  for  both  cases.  A  severe  log-law  mismatch  is  not  observed  in 
any  case,  but  the  velocity  profile  is  not  completely  smooth,  exhibiting  a  small  hump  that  is 
reduced  in  size  as  the  time-averaged  transition  position  shifts  further  away  from  the  wall.  A 
lower  value  of  Cs  also  seems  to  promote  a  smoother  profile.  In  all  cases,  the  computed  profiles 

in  the  logarithmic  and  wake  regions  lie  above  the  theoretical  solution.  This  is  a  consequence  of  a 
slight  under-prediction  of  the  wall  shear  stress  exhibited  by  all  models  tested. 
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a.)  Luker,  et  al.  b.)  Smits  and  Muck 

Figure  3.3.1:  Effect  of  ensemble-averaging  technique  and  Cs 

3.4.  Comparison  with  IDDES 


An  alternative  wall-modeled  LES  strategy  (termed  IDDES)  has  been  developed  by  Shur,  et  al. 
[10]  as  an  extension  of  the  delayed  detached-eddy  simulation  (DDES)  model.  [6]  Variants  based 
on  the  Menter  SST  model  and  the  Spalart-Allmaras  model  have  been  formulated.  The  general 
idea  adopted  by  IDDES  is  that  the  mesh  spacing,  plus  the  presence  of  initial  fluctuation  content, 
will  dictate  whether  the  model  responds  as  a  wall-modeled  LES,  as  DDES,  or  as  a  pure  RANS 
model.  It  is  thus  a  more  general  framework  than  our  model,  which  is  strictly  a  type  of  wall- 
modeled  LES.  The  wall-modeled  LES  branch  of  IDDES  differs  from  ours  in  several  respects. 
First,  the  Menter-SST  model  is  used  as  the  basis,  rather  than  the  Menter  BSL  model.  Secondly, 
only  the  destruction  term  in  the  turbulence  kinetic  energy  equation  is  modified  through  the  action 


y/5  (bottom);  log10(y/5)  (top)  y/5  (bottom);  log10(y/S)  (top) 


Figure  3.4.1:  Reynolds-shear  stress,  turbulence 
kinetic  energy,  and  mean  velocity  predictions 
for  different  turbulence  models  (Elena  and 
Lacharme  experiment) 


Figure  3.4.2:  Reynolds-shear  stress,  turbulence 
kinetic  energy,  and  mean  velocity  predictions 
for  different  turbulence  models  (  Luker,  et  al. 
experiment) 
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of  a  blending  function.  The  eddy  viscosity  responds  to  the  combined  influences  of  the 
turbulence  kinetic  energy,  the  turbulence  frequency,  and  a  specified  filter  width  -  there  is  no 
separate  SGS  eddy  viscosity.  In  our  model,  the  eddy  viscosity  description  is  designed  to  shift 
between  RANS  and  LES  branches,  and  in  principle,  different  SGS  and  RANS  models  could  be 
employed.  The  blending  functions  used  in  IDDES  are  fundamentally  based  on  the  ratio  of  the 
wall  distance  to  a  filter  width,  and  thus  the  transition  between  RANS  and  LES  branches  occurs  at 
fixed  locations  for  a  particular  grid.  In  our  model,  the  instantaneous  transition  location  will 
migrate  in  response  to  the  resolved  and  ensemble-averaged  solution.  The  IDDES  method  uses 
flow-dependent  ‘enhancement’  functions  to  augment  the  RANS  component  in  the  vicinity  of  the 
transition  location  -  the  authors  argue  that  such  functions  are  crucial  for  combating  Tog-law 
mismatch’.  For  the  array  of  flat-plate  boundary  layers  considered  in  this  work,  the  transition 
location  as  calculated  in  the  IDDES  model  is  at  y/8o=l/40.  In  wall  coordinates,  this  places  the 
transition  at  d+  =  17.8,  36.7,  and  251.8  for  the  Elena  and  Lacharme,  Luker,  et  al,  and  Smits  and 
Muck  cases,  respectively  The  first  two  locations  are  near  the  beginning  of  the  logarithmic 
region,  while  the  latter  is  more  toward  the  middle  of  the  region. 

Figures  3.4.1-3.4.3  compare  solutions  obtained  using  the  new  model  (Cn  =  1.5,  Cs  =  15)  with 
those  from  IDDES  and  with  Menter  BSL  RANS  solutions  for  each  of  the  three  cases.  The 
RANS  velocity  solutions  agree  very  well  with  the  compressible  law  of  the  wall  and  with  the 
experimental  data  and  thus  might  be  used  as  a  reference  to  assess  the  predictive  capability  of  the 
other  models.  The  resolved  turbulence  kinetic  energy  and  Reynolds  shear-stress  are  plotted  for 
the  LES/RANS  models,  while  the  modeled  turbulence  kinetic  energy  and  Reynolds  shear-stress 
are  plotted  for  the  RANS  model.  The  axial  velocity  is  plotted  versus  logio(y/8o)  and  is  scaled  by 
the  free-stream  velocity  instead  of  the  friction  velocity.  This  removes  ambiguities  associated 
with  the  possible  under-prediction  of  the  wall  shear  stress.  Both  LES/RANS  models  predict 
outer-layer  values  of  turbulence  kinetic  energy  and  Reynolds-shear  stress  that  are  comparable  to 
those  modeled  using  RANS.  The  model  responses  differ  greatly  in  the  inner  layer,  with  the 
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IDDES  solution  exhibiting  a  peak  in  turbulence  kinetic  energy  also  captured  in  the  CN=  1.2 

solution  of  Figure  3.2.1.  The  rapid  diminishing  of  the  resolved  fluctuation  energy  near  the  wall 
by  the  new  model  is  again  in  evidence.  The  new  LES/RANS  model  predicts  velocity  profiles 
that  are  in  close  accord  with  the  RANS  solutions,  whereas  the  IDDES  solutions  show  an  increase 
in  velocity  toward  the  outer  part  of  the  logarithmic  region  followed  by  a  decrease  in  velocity  in 
the  inner  part  of  the  layer.  Figure  3.4.4  presents  a  scatter  plot  of  predicted  skin  friction  versus 
CN  for  all  of  the  models  considered  in  this  study.  For  IDDES,  CN  is  arbitrarily  set  to  one, 

while  for  Choi,  et  al,  CN  is  set  to  1.5.  The  general  trend  is  that  all  of  the  LES/RANS  models 

under-predict  the  wall  shear  stress 
level  (as  provided  by  the  RANS 
model),  with  the  IDDES  and  CN  = 

1.2  solutions  showing  the  largest 
deviations.  Based  on  these  results 
and  the  others  presented,  it  appears 
that  a  RANS-to-LES  transition  that 
occurs  near  the  inner  edge  of  the 
logarithmic  layer  can  result  in 
insufficient  dissipation  of  resolved 
turbulence  energy  in  the  outer  part 
of  the  logarithmic  layer.  Axial 
fluctuation-intensity  predictions 
improve,  but  inaccuracies  in  the 
wall  shear  stress  and  inner-layer 
velocity  profile  appear.  Shifting 
the  RANS-to-LES  transition 
toward  the  middle  and  outer  parts 
of  the  logarithmic  layer,  as  is  done 
in  the  new  model  with  CN  =1.5-1. 6, 

yields  closer  agreement  with  theory  and  with  RANS  predictions  for  the  velocity  profile  but  leads 
to  an  under-prediction  of  the  resolved  axial  fluctuation  intensity  in  the  near- wall  region.  Figure 
3.4.5  shows  that  the  differences  between  the  predictions  of  the  new  model  and  those  of  IDDES 
are  not  likely  to  be  a  result  of  the  SGS  component  of  the  closure.  The  normalized  eddy  viscosity 
is  much  larger  for  IDDES  in  the  outer  layer,  where  the  SGS  model  should  be  most  active.  The 
opposite  situation  occurs  close  to  the  wall.  The  degree  to  which  these  observations  are  also  a 
function  of  the  behavior  of  the  current  numerical  method  is  unknown,  and  the  sensitivity  of  the 
predictions  to  mesh  refinement  remains  to  be  assessed.  It  is  unlikely,  however,  that  practical 
LES/RANS  calculations  using  meshes  finer  (in  terms  of  cell  size  per  shear  layer  thickness)  than 
those  considered  herein  will  be  affordable  in  the  near  term,  but  the  methods  should  be  assessed 
for  coarser  levels  of  resolution  (10  to  15  cells  /  shear  layer  thickness). 


y/50 

Figure  3.4.5:  Normalized  kinematic  eddy  viscosity 
distributions  -  new  model  and  IDDES 
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4.  Model  Assessment:  16  degree  Smooth  Compression  Ramp 

The  LES/RANS  models  have  been  applied  to  a  Mach  2.86,  16-degree  smooth  compression- 
comer  interaction  mapped  by  Donovan,  et  al.  [1 1,12]  The  flow  conditions  for  this  case  and  the 
previously-discussed  flat-plate  cases  are  shown  in  Table  4.1.  The  specific  variants  tested  are  the 
original  model  of  [2]  (with  the  model  constant  al  =  13.06  5 +  76.7  fitted  as  a  function  of  the 
surface  coordinate  s ),  the  new  model  with  CN  =1.5  and  Cs  =  1 5  using  Weighing  #1  for 
ensemble  averaging,  and  the  new  model  with  CN  =1.5  and  Cs  =  1 5  using  Weighing  #2  for 

ensemble  averaging  (See  Section  2.4).  The  mesh  extends  from  X  =  -0.28  m  to  X  =  0.2717  m  in 
the  streamwise  direction,  with  X  =  0  m  corresponding  to  the  start  of  the  16  degree  turn.  The 
mesh  extends  0.14  m  in  the  wall-normal  (Y)  direction  and  +/-  0.0756  m  in  the  spanwise  (Z) 
direction.  The  mesh  resolution  is  such  that  20  cells  /  incoming  boundary  layer  thickness  are 
present  in  the  wall-transverse  directions,  and  the  total  number  of  interior  mesh  cells  is  8.64  M. 
The  recycle  plane  is  located  7.5  boundary-layer  thicknesses  (S=  2.5  cm)  downstream  of  the 
inflow  plane. 


Table  4.1:  Free-stream  and  boundary-layer  properties  for  experiments  considered  in  this  study. 


Case 

S0 ,  mm 

Uo o ,  m/s 

Res 

Pa  >  Pa 

T,  K 

cf 

(xlO  3) 

Elena  and  Lacharme  [21] 

2.32 

10  (12)* 

552 

5.59xl04 

5.0xl04 

291 

2.15 

Luker,  et  al.  [22] 

2.80 

9.9 

602 

1.78xl05 

2.1xl05 

298 

1.6 

Smits  and  Muck  [23] 

2.79 

25 

562 

1.58xl06 

6.9xl06 

263 

1.07 

Donovan  ,  et  al.  [16] 

2.86 

28 

580 

1.76xl06 

6.9xl06 

270 

1.08 

*-  evaluated  at  99.9%  of  boundary  layer  edge  velocity 


Figure  4.1.1:  Iso-surfaces  of  swirl  strength  (7500  s-1)  colored  by  temperature  values  (Donovan,  et  al. 
experiment) 
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Figure  4.1.2:  Surface  wall  pressure  distributions 
(Donovan,  et  al.  experiment) 
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Figure  4.1.3:  Surface  skin  friction 
distributions  (Donovan,  et  al.  experiment) 


An  iso-surface  of  swirl  strength  (7500  s'1),  colored  by  temperature,  is  shown  in  Figure  4.1.1  for 
the  Donovan,  et  al.  experiment.  The 


swirl  strength  is  defined  as  the 
imaginary  portion  of  the  complex 
eigenvalue  of  the  local  velocity 
gradient  tensor.  The  higher  the  value 
of  the  swirl  strength,  the  shorter  the 
time  required  for  a  fluid  particle  to 
swirl  about  a  vortex  core.  Higher 
values  of  the  swirl  strength  typically 
correspond  to  smaller-scale  turbulent 
structures,  and  the  distribution  of  the 
swirl  strength  can  provide  an 
approximate  measure  of  the 
distribution  of  the  sizes  of  turbulent 
eddies.  The  value  chosen  (7500  s'1) 
is  large  enough  to  highlight  the 
presence  of  longitudinally-oriented 
vortical  structures  in  the  recovering 
boundary  layer  downstream  of  the 
isotropic  compression.  These  appear 
to  arise  from  amplification  of 
inhomogeneities  in  the  incoming 
boundary  layer  and  have  been 
observed  in  both  experiments  and 
computations  of  compression-ramp 
flows. 


V/5ref 

Predictions  for  mean- flow  properties  Figure  4.1.4:  Streamwise  mass  flux  distributions 

for  the  Donovan,  et  al.  experiment  throughout  the  interaction  region  (Donovan,  et  al. 

are  shown  in  Figures  4. 1.2-4. 1.4.  experiment) 
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Time-  and  span-averaging  of  the  instantaneous  data  is  used  to  generate  the  mean  values.  Surface 
pressure  distributions  in  Figure  4.1.2  show  evidence  of  a  smooth  compression  process  unaffected 
by  local  flow  separation.  Figure  4.1.3  shows  that,  relative  to  the  original  model,  the  new 
LES/RANS  model  yields  better  agreement  with  surface  skin  friction  measurements  of  Donovan 
within  and  downstream  of  the  compression  region. 

Figure  4.1.4  compares  mean  mass  flux  distributions  (based  on  the  velocity  component  tangential 
to  the  surface)  with  experimental  hot-wire  data  throughout  the  interaction.  Good  agreement  is 
indicated  in  general,  though  the  computed  profdes  do  show  an  under-prediction  of  the  mass  flux 
distributions  in  the  near-wall  region  downstream  of  the  compression  region.  Locations  in  the 
flow  field  associated  with  the  coalescence  of  compression  waves  into  a  shock  wave  (Figure  19, 
X  =  0.1016  and  X  =  0.1278  stations)  are  also  not  well-predicted.  The  hot-wire  measurements 
may  not  be  as  accurate  in  this  region.  The  new  model  predicts  a  faster  recovery  of  the  inner  part 


y/s0 

Figure  4.1.5:  Mass  flux  fluctuation  intensity 
profiles  at  different  streamwise  stations 
(Donovan,  et  al.  experiment) 
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Figure  4. 1 .6:  Reynolds  axial  stress  profiles  at 
different  streamwise  stations  (Donovan,  et  al 
experiments) 
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of  the  boundary  layer  downstream  of  the  compression  region  and  provides  slightly  better 
agreement  with  experiment  in  this  region. 

Figure  4.1.5  compares  mass-flux  fluctuation  intensity  predictions  with  experimental  data  at  three 
streamwise  locations  within  the  interaction  region.  The  mass-flux  fluctuation  intensity  is 
computed  according  to 

^,)']2)1/2  =[(/^)2-^,^]1/2  (4.1.1) 

where  the  overbar  represents  time  and  span-averaging  of  the  grid-filtered  data  and  wyis  the 

velocity  vector  expressed  in  a  coordinate  system  that  is  aligned  with  the  compression-comer 
surface.  At  its  peak  in  the  incoming  boundary  layer,  the  maximum  value  of  the  mass  flux 
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Figure  4.1.7:  Reynolds  normal  stress  profiles 
at  different  streamwise  stations  (Donovan,  et 
al.  experiments) 
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Figure  4.1.8:  Reynolds  shear  stress  profiles  at 
different  streamwise  stations  (Donovan,  et  al. 
experiments) 


fluctuation  intensity  is  about  15%,  increasing  to  a  maximum  of  about  25%  midway  along  the 
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compression-ramp  portion.  As  the  boundary  layer  begins  to  relax  back  to  a  new  equilibrium 
state,  the  mass  flux  fluctuation  intensity  diminishes.  The  LES/RANS  predictions  are  in  good 
agreement  with  experimental  data  at  every  station,  with  the  new  model  providing  better  overall 
agreement  nearer  to  the  surface. 


Favre-averaged  Reynolds-stresses  are  computed  according  to 


— ~  rr^  tt  ~  ~ 

piijU  =  putu  - 


pu,Plij 

P 


(4.1.2) 


Only  the  resolved  Reynolds  stresses  are  captured  in  this  approach,  but  as  the  experimental  data 
does  not  extend  far  into  the  wall  layer,  this  is  adequate  for  comparison.  Figure  4.1.6  compares 
experimental  measurements  of  the  Reynolds  axial  stress  with  computational  predictions.  Good 
agreement  is  generally  indicated,  though  calculations  performed  using  the  original  LES/RANS 
model  do  over-predict  the  peak  value  of  the  Reynolds  normal  stress  at  the  X  =  0.1516  m  data 
station,  which  is  near  the  point  where  the  isentropic  compression  ends.  The  new  model  provides 
improved  results  at  this  station.  Evidence  obtained  from  earlier  shock  impingement  and 
compression-ramp  studies  indicates  that  the  original  model  may  over-amplify  turbulent 
fluctuations  in  the  near-wall  region  downstream  of  a  strong  compression.  This  problem  may  be 
related  to  a  collapse  in  the  blending  function  toward  the  wall,  which  enables  larger  turbulent 
eddies  to  interact  with  the  wall  without  significant  attenuation.  The  fact  that  the  new  model 
improves  upon  this  response  is  encouraging. 


Comparisons  with  the  measured  Reynolds  normal  stress  (Figure  4.1.7)  show  that  the  trends 
regarding  amplification  of  this  component  are  in  agreement  with  those  shown  in  the 
experimental  data  and  are  consistent  with  those  evidenced  in  Figures  4.1.5  and  4.1.6. 
Agreement  with  experiment  is  poor  for  the  Reynolds  shear  stress  (Figure  4.1.8),  as  the 
computations  over-predict  the  measured  values  by  more  than  a  factor  of  two  at  the  most 
downstream  stations.  The  experimental  data  shows  a  marked  decrease  in  the  amplification  rate 
of  the  Reynolds  shear  stress  downstream  of  the  end  of  the  compression  region.  The  calculations 
do  not  show  this  response,  and  neither  the  experimental  data  nor  the  computational  predictions 
for  the  other  Reynolds-stress  components  show  such  a  rapid  return  to  equilibrium. 


5.  Model  Assessment:  20  degree  Sharp  Compression  Ramp 

5.1.  Mean  flow  predictions 

This  section  describes  the  application  of  the  new  LES/RANS  model  to  Mach  3  flow  over  a  20- 
degree  compression  comer  also  studied  at  Princeton  University.  Velocity,  wall-pressure,  and 
skin  friction  measurements,  along  with  hot-wire  measurements  of  rms  mass  flux,  Reynolds  axial 
stress,  and  Reynolds  shear  stress  are  available  for  this  case  in  the  Supersonic  /  Hypersonic  Shock 
Wave  /  Turbulent  Boundary  Layer  Interaction  Database  [13].  This  case  differs  from  the 
Donovan  and  Smits  experiment  in  that  the  oblique  shock  generated  at  the  compression  comer  is 
strong  enough  to  separate  the  incoming  turbulent  boundary  layer.  We  have  applied  the  new 
model  (CAf=  1.5,  Cs  =15)  to  this  case  using  ensemble-averaging  weighting  methods  1  and  2  as 

well  as  the  original  approach  of  Choi,  et  al.  Figure  5.1.1  shows  snapshots  of  temperature 
contours  for  the  Choi,  et  al.  model  (PPM  discretization)  and  the  new  model  (LD-PPM,  weighing 
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method  2).  Time-averaged  wall-pressure  and  skin  friction  distributions  are  shown  in  Figure 
5.1.2.  The  new  models  (‘newest  model’  corresponds  to  the  use  of  weighting  method  1)  give 
essentially  identical  results  and  provide  a  slight  improvement  over  the  Choi,  et  al.  model  in  the 
prediction  of  the  upstream  extent  of  axial  separation.  Still,  all  models  appear  to  over-predict  the 
extent  of  separated  flow  -  a  fact  further  evidenced  in  the  axial  velocity  profiles  of  Figure  5.1.3. 
The  wake-like  response  encountered  as  the  boundary  layer  moves  over  the  low-momentum  flow 
region  at  the  corner  is  over-predicted  by  all  models.  This  leads  to  a  delayed  recovery  of  the 
turbulent  boundary  layer  downstream  of  re-attachment. 
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Figure  5.1.1:  Centerplane  temperature  contours  (top:  LD-PPM,  new  model;  bottom: 
PPM,  model  of  Choi,  et  al. 


Figure  5.1.2:  Surface  pressure  and  skin  friction  distributions:  20-degree  compression  comer  24 


-  Choi  et.  al.  (2009) 

- -  NewModel 


-  NewModel  (LD-PPM) 

-  Newest  Model 


5.2.  Reynolds  stress  predictions 


Predictions  of  the  rms  mass  flux  fluctuation  intensity  and  Reynolds  stress  fields,  with  the  latter 
estimated  both  by  Favre  ensemble-averaging  and  by  the  use  of  Morkovin’s  strong  Reynolds 
analogy  (SRA)  as  utilized  by  Smits  and  Muck,  are  described  in  this  section.  Considering  only 
the  Reynolds  axial  and  shear  stresses,  Favre-averaging  yields 


and 


pu[u[  =  puxux  - 


puxpux 

P 


(Reynolds  axial  stress ) 


(5.2.1) 


— j  -  pUxpU2 

puxu2  =  puxu 2 - = —  (Reynolds  shear  stress) 

P 

Morkovin’s  SRA  requires  first  an  expression  for  the  rms  mass  flux  fluctuation  intensity: 


(5.2.2) 


{(pu,)']2)'2  =[(pux)2 -puxpux]U2  (5.2.3) 

from  which 

pu[u[  =  — ^ -  (Reynolds  axial  stress)  (5.2.4) 

P P 
and 

pu[u'2  =  ^( pii\)'u'2  (Reynolds  shear  stress)  (5.2.5) 

can  be  obtained.  Here,  p2  =  \-2RuT(y -\)MX2  +  [(j -\)M2]2 ,  with Mx being  the  local 
average  Mach  number  evaluated  as  ux  /  c ,  where  c  is  the  local  sound  speed,  and  RuT  is  the 

correlation  coefficient  u\T' /\(u\)2 1  U7’')2 )  ,  which  is  set  to  -0.8.  Figure  5.2.1  shows  the 
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evolution  of  the  rms  mass  fluctuation  intensity  throughout  the  interaction,  while  Figure  5.2.2 
shows  the  evolution  of  the  Reynolds  axial  and  shear  stress  distributions.  The  Reynolds  axial 
stress  is  most  in  error  in  the  low-momentum  region  located  near  the  wall  downstream  of  re¬ 
attachment.  Morkovin’s  SRA  provides  predictions  that  are  close  to  those  provided  by  Favre- 
averaging,  and  the  fact  that  the  rms  fluctuation  intensity  is  actually  over-predicted  by  the 


Figure  5.2.1:  Mass  flux  fluctuation  intensity  distributions 

LES/RANS  methods  means  that  the  most  probable  cause  of  the  discrepancy  is  an  over-prediction 
of  the  vertical  extent  of  the  low-momentum  region,  as  the  average  density  and  /? 2  would  both  be 
lower  than  in  the  experiment.  The  Reynolds-shear  stress  predictions  are  off  by  at  least  a  factor 
of  two,  but  the  uncertainty  in  the  experimental  measurements  due  to  probe  misalignment  and 
calibration  errors  is  believed  to  be  very  high. 

5.3.  Shock  motion:  20  degree  interaction 

We  have  used  several  sampling  techniques  to  study  aspects  of  separation-shock  motion  for  the 
20-degree  interaction.  There  appear  to  be  two  distinct  responses  present,  both  of  which  can 
influence  the  observed  intermittency  characteristics  in  this  flow.  As  shown  in  Figure  5.3.1, 
streak-like  structures  present  in  the  incoming  boundary  layer  can  disturb  the  shock  front,  pushing 
it  upstream  and  downstream  of  its  time-averaged  position.  The  extent  of  shock  motion  can  be 
quantified  by  an  intermittency  distribution,  defined  as  the  amount  of  time  that  the  shock  is 
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Smits  and  Muck  [5] 
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New  model  (CN  =  1.5,  Weighting  #1) 


Smits  and  Muck  [5] 
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Figure  5.2.2:  Reynolds  axial  stress  (left)  and  shear  stress  (right)  distributions 


upstream  of  a  particular  location.  This  is  shown  in  Figure  5.3.2  for  the  20-degree  interaction. 
The  shock  front  moves  over  a  range  of  -0.6  boundary  layer  thicknesses.  A  low-frequency  signal 


Y/6=0.2  Y/6=0.7 


Figure  5.3.1:  Velocity  magnitude  contours  at 
different  wall-normal  planes 


Figure  5.3.2:  Intermittency  and  probability  of 
re  versed-flow  distributions. 
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is  also  observed,  as  indicated  in  the  spectral  maps  shown  in  Figure  5.3.3.  Here,  the  normalized 
power  spectrum  is  plotted  at  different  X  stations  throughout  the  interaction.  A  shift  to  a  lower- 
frequency  band  is  noted  near  the  separation  shock  position,  followed  by  a  recovery  toward  a 
typical  boundary  layer  spectrum  as  the  flow  re-attaches  on  the  compression  ramp.  Also  noted  in 
Figure  5.3.3  are  two  estimates  for  the  dominant  low-frequency  signal,  one  from  the  correlation  of 
Dupont  [14]: 

Sr  =  ^'Lsep  ,Sr  =0.025,  (5.3.1) 

^00 

and  another  from  a  residence-time  distribution  analysis  to  be  discussed  later.  Both  estimates  are 
in  close  agreement  with  the  observed  peak  in  the  low-frequency  signal. 
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Figure  5.3.3:  Spectral  maps  extracted  at  different  Y 
locations  -  20-degree  interaction. 

6.  Model  Assessment:  24-degree  Sharp  Compression  Ramp 

6.1.  Mean  flow  behavior 

Here,  we  repeat  many  of  the  analyses  performed  for  the  20-degree  interaction  for  the  case  of 
nominal  Mach  3  flow  over  a  24  degree  compression  comer.  [13]  The  shock  wave  in  this  case  is 
strong  enough  to  induce  a  large  region  of  separated  flow  -  extending  ~6  cm  upstream  of  the 
wedge  apex.  Figures  6.1.1  and  6.1.2  present  surface  pressure  and  skin  friction  distributions, 
respectively.  As  noted,  all  LES/RANS  models  over-predict  the  upstream  extent  of  flow 
separation  for  this  case.  The  newer  models,  perform  better  in  this  regard,  and  there  is  little 
difference  in  the  predictions  with  respect  to  the  choice  of  the  ensemble-averaging  technique. 
Figure  6.1.3  shows  velocity  profiles  extracted  throughout  the  interaction.  These  also  indicate  an 
over-prediction  of  the  size  of  the  backflow  region  and  an  associated  delay  in  the  boundary-layer 
recovery  rate. 
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Figure  6.1.1:  Surface  pressure  distributions  Figure  6.1.2:  Surface  skin  friction  distributions 


Choiet.  al.  (2009) 

NewModel  (C  H  =  1  .5,  Weighting  #1 ) 
New  Mo  del  (C  N  =  1  .5,  Weighting  #2) 


6.2.  Reynolds  stress  predictions. 

As  in  the  20-degree  interaction,  we  have  extracted  Reynolds  stress  statistics  using  Favre 
ensemble  averaging  and  using  Morkovin’s  strong  Reynolds  analogy.  Predictions  of  rms  mass 
flux  fluctuation  intensity  in  Figure  6.2.1  show  good  agreement  with  the  peak  values  and  their 
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(strong  Reynolds  analogy) 

locations,  but  this  is  somewhat  misleading,  as  the  intensity  curves  are  normalized  with  respect  to 
the  mean  momentum  flux,  which  also  is  in  error  due  to  the  separation  prediction  being  incorrect. 
As  such,  the  Reynolds  axial  stress  distributions  (Figure  6.2.2)  show  significant  discrepancies, 
and  the  use  of  Morkovin’s  SRA  provides  no  real  improvement.  These  results  provide  more 
evidence  that  inaccuracies  in  the  prediction  of  the  structure  of  the  backflow  region  affect  the 
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Figure  6.3.1:  Power  spectra  distributions:  24-degree  ramp 
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entirety  of  the  downstream  flow. 


6.3.  Shock  motion  -  24  degree  interaction  and  residence  time  analysis 


We  have  also  extracted  wall  pressure  spectra  for  the  24  deg.  case,  and  as  expected,  it  also  shows 
a  shift  to  lower  frequencies  in  the  vicinity  of  the  separation  shock.  Figure  6.3.1  shows  wall 
pressure  spectra  for  the  24-degree  ramp.  A  larger  region  of  low-frequency  motion  is  indicated 
upstream  of  the  compression  ramp,  indicating  that  the  response  scales  with  the  size  of  the 
separation  region. 


Figure  6.3.2:  Streamlines  entering  separation  Figure  6.3.3:  Streamlines  entering  separation 

region  (view  from  above,  24-degree  interaction )  region  (view  from  above,  20-degree  interaction) 


The  lines  corresponding  to  “Residence  time  analysis”  and  “Residence  Time  Frequency”  in  the 
above  figures  are  from  a  new  analysis  that  extracts  residence-time  distributions  for  fluid  particles 
entering  the  recirculation  zones.  The  fact  that  a  large  amount  of  experimental  data  can  be 
collapsed  into  a  relatively  narrow  range  of  Strouhal  numbers  according  to  Dupont’s  correlation 
[14]  suggests  that  the  dominant  time  scale  must  be  one  associated  with  entrainment  of  fluid  into 
and  out  of  the  recirculation  region,  which  naturally  would  scale  with  the  length  of  the  separation 
region.  To  explore  this,  a  collection  of  4000  evenly-distributed  streamlines  passing  through  the 
separation  region  were  analyzed  to  determine  the  amount  of  time  a  fluid  particle  would  spend  in 
this  region.  To  do  this,  a  line  integral  of  velocity  magnitude  was  performed  along  each 
streamline,  starting  when  the  streamline  crossed  x  =  -5  cm  and  ending  once  it  had  crossed  x  =  2 
cm.  Figure  6.3.2  shows  a  top  view  of  an  ensemble  of  near-surface  streamlines  for  the  24-degree 
interactions,  while  Figure  6.3.3  shows  the  same  for  the  20-degree  interaction.  A  weakly- 
separated  flow,  characterized  by  significant  spanwise  migration  of  fluid  particles  and  complex 
topological  features,  is  found  in  the  region  of  intermittent  separation- shock  motion.  Further 
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downstream,  a  more  ‘two-dimensional’  separation  pattern  is  found  at  intervals  in  the  spanwise 
(Z)  direction.  The  blue  lines  in  Figure  6.3.2  indicate  the  X  stations  where  integration  begins  and 
ends.  The  red  rectangle  highlights  streamlines  within  one  of  the  ‘two-dimensional’  sections  of 
the  separated  region,  where  the  vorticity  vector  points  predominately  in  the  Z  direction. 

Once  the  times  were  calculated  for  each  streamline,  a  probability  density  function  for  the 
residence  time  distribution  was  created,  as  shown  in  Figure  6.3.4.  The  most  probable  residence 
time  is  3.23e-3  seconds  for  the  24-degree  interaction  and  2.25e-3  for  the  20-degree  interaction. 
The  frequencies  associated  with  these  residence  times  are  close  to  the  shock  motion  frequencies 
predicted  for  this  case  by  the  Dupont  et  al.  correlation  and  coincide  with  one  of  the  local  peaks  in 
the  spectral  map.  These  results  provide  evidence  of  a  connection  between  the  observed  low- 
frequency  signal  and  the  time  required  for  fluid  to  enter  and  leave  the  recirculation  zone.  The 
fact  that  residence  time  distributions  of  this  type  can  be  obtained  from  mean-flow  data  indicates 
that  it  might  be  possible  to  predict  a  dominant  low-frequency  model  without  conducting  an 
unsteady  analysis.  The  challenge  is  in  identifying  entrainment  pathways  so  that  an  appropriate 
sampling  window  can  be  constructed. 


Figure  6.3.4:  Residence  time  distributions  (20-degree  and  24-degree 
ramps 


7.  Model  Assessment:  3D  Shock  /  Boundary  Layer  Interaction 

We  have  conducted  LES/RANS  simulations  of  a  Mach  2.5  shock  /  boundary  layer  interaction  in 
a  wind  tunnel  (experiments  conducted  at  Cambridge  University  [15])  as  a  means  of  assessing 
methods  described  in  Section  2.3  for  maintaining  a  RANS-type  response  in  regions  of  the  flow 
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that  are  not  resolved  enough  to  sustain  turbulence.  The  shock  is  generated  by  an  8-degree 
compression  ramp  placed  on  the  top  of  the  wind  tunnel.  The  shock  impinges  upon  the  bottom 
surface  of  the  wind  tunnel,  creating  a  region  of  shock-separated  flow.  The  structure  of  the  SBLI 
is  known  to  be  impacted  by  comer  vortices  that  are  generated  through  interactions  between  the 
oblique  shock  and  the  sidewall  boundary  layers.  These  force  fluid  toward  the  comer,  reducing 
the  effective  cross-sectional  area  of  the  wind  tunnel  and  causing  the  core  flow  to  accelerate 
around  the  displaced  viscous  layer.  The  mesh  contains  about  46  M  cells  but  is  not  refined 
significantly  in  the  Y  direction,  rendering  it  unsuitable  for  sustaining  turbulence  in  the  sidewall 
boundary  layers.  Figure  7.1.1  shows  a  center-plane  snapshot  of  temperature,  illustrating  the  fact 
that  turbulence  is  sustained  in  the  top  and  bottom  surface  boundary  layers.  Thickening  of  the 
boundary  layer  due  to  shock  impingement  is  also  indicated.  The  impact  of  comer  vortices  on 
the  near-surface  flow  field  is  illustrated  in  Figure  7.1.2,  a  plot  of  near-surface  streamlines  colored 
by  axial  velocity  magnitude.  The  left  figure  corresponds  to  the  use  of  the  ‘RANS  wall’  option, 
while  the  right  figure  corresponds  to  the  use  of  grid-scale  limiting  (see  Section  2.3)  The 
predictions  are  similar,  indicating  that  the  structure  of  the  comer  vortices  is  not  significantly 
impacted  by  this  choice. 
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Figure  7.1.1:  Centerplane  temperature  contours:  Cambridge  SBLI 

The  effect  of  both  techniques  in  attenuating  resolved-scale  turbulence  is  illustrated  in  Figure 
7.1.3,  which  plots  instantaneous  temperature  (bottom)  and  eddy  viscosity  (top)  contours  at  the 
shock  crossing  location.  The  grid-scale  limiting  method  is  somewhat  less  aggressive,  as  it 
produces  smaller  eddy  viscosities  and  does  allow  some  turbulence  content  to  remain  in  regions 
near  the  comers. 

Centerline  velocity  profiles  within  the  interaction  and  centerline  surface  pressure  distributions 
are  shown  in  Figure  7.1.4.  These  indicate  that  the  LES/RANS  method  performs  very  well  in 
capturing  the  response  of  the  viscous  layer  to  the  shock  wave. 
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X  X 

Figure  7.1.2:  Near-wall  flow  structure  for  Cambridge  SBLI  (left:  ‘RANS  wall’;  right:  grid-scale 
limiting) 


Z  (m)  Z  (m) 

Figure  7.1.3:  Eddy  viscosity  (top)  and  temperature  (bottom)  contours  at  shock-crossing  location  (left: 

‘RANS  wall’;  right:  grid-scale  limiting) 
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Figure  7. 1 .4:  Centerline  velocity  profiles  (left)  and  surface  pressure  distribution  (right):  Cambridge 
SBLI 


8.  Model  Assessment:  Airfoil  near  Static  Stall 


It  is  of  interest  to  investigate  the  performance  of  the  LES/RANS  model  for  low-speed  flows 
characteristic  of  helicopter  aerodynamics,  as  there  is  a  significant  interest  within  ARO  in  better 
prediction  methods  for  dynamic  stall.  In  our  past  work,  we  have  used  recycling  /  rescaling 
techniques  to  sustain  turbulence  within  boundary  layers.  These  will  not  generally  work  for 
airfoil-type  flows,  as  they  rely  on  a  RANS  base  state  and  as  airfoil  turbulent  boundary  layers  will 
be  initiated  by  growth  of  natural  instabilities  within  a  laminar  boundary  layer  or  within  a  laminar 


Figure  8.1.1:  X-Y  centerplane  mesh  for  A- 
airfoil 


Figure  8.1.2:  Iso-surfaces  of  swirl  strength 
(2000  s'1)  illustrating  development  of  eddy 
structures  in  airfoil  boundary  layer 


separation  region.  It  is  of  interest,  therefore,  to  determine  the  baseline  response  of  the  current 
LES/RANS  methods  for  flow  over  an  airfoil.  To  this  end,  we  have  conducted  a  simulation  of 
flow  over  an  ‘A- Airfoil’  (an  Aerospatiale  design)  at  conditions  near  static  stall  [16].  Mary  and 
Sagaut  [17]  previously  conducted  a  large-eddy  simulation  of  this  flow.  The  free-stream  Mach 
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number  is  0.15,  the  Reynolds  number  based  on  a  chord  length  of  0.6  m  is  2.1e6,  and  the  angle  of 
attack  is  13.3  degrees.  An  X-Y  snapshot  of  the  computational  mesh,  which  contains  ~30e6  cells, 


Figure  8.1.3:  Streamwise  velocity  profiles  along  airfoil  surface. 


is  shown  in  Figure  8.1.1,  and  an  iso-surface  of  swirl-strength,  illustrating  the  growth  of  large 
turbulent  eddies  on  the  suction  side,  is  shown  in  Figure  8.1.2.  Comparisons  with  experimental 
velocity  profile  data  (Figure  8.1.3)  and  turbulence  intensity  data  (Figure  8.1.4)  are  generally 
favorable,  though  the  turbulence  intensity  level  is  under-predicted  near  the  trailing  edge,  an 
effect  also  observed  by  Mary  and  Sagaut  [16].  The  first  model  tested  (Choi,  et  al.  [2)  requires  a 


Figure  8. 1 .4:  Streamwise  rms  velocity  fluctuation  profiles  along  airfoil  surface 


pre-selection  of  a  model  constant  that  varies  with  chord  along  the  upper  and  lower  surfaces.  The 
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calculation  of  this  quantity,  described  in  detail  in  [2],  requires  an  estimate  of  the  boundary  layer 
edge  and  the  flow  properties  at  that  edge.  This  information  was  obtained  for  the  lower  and  upper 
airfoil  surfaces  from  an  initial  RANS  calculation.  The  new  model,  developed  during  this  grant, 
does  not  need  this  information.  Comparisons  with  surface  pressure  coefficient  (Figure  8.1.5)  and 
skin  friction  (Figure  8.1.6)  distributions  also  show  good  agreement  with  experiment.  However, 
the  reference  for  the  experiment  [15]  and  the  LES  computations  presented  in  [16]  indicate  that  a 
laminar  separation  bubble  should  be  present.  Neither  of  the  LES/RANS  models  predicts  a 
region  of  laminar  flow  near  the  leading  edge,  indicating  the  need  to  include  a  transition  model 
for  the  RANS  component  of  this  closure. 


Figure  8.1.5:  Surface  pressure  coefficient  for  Figure  8.1 .6:  Surface  skin  friction  coefficient  for 

Aerospatiale  ‘A-airfoil’  Aerospatiale  ‘A-airfoil’ 


9.  Large-eddy  Simulations  of  the  Elena-Lacharme  Flat-Plate  Boundary  Layer 
Experiment 

Though  the  LES/RANS  model  works  reasonably  well,  there  are  still  aspects  of  the  formulation 
that  need  improvement.  Among  these  is  the  fact  that  accurate  predictions  of  the  mean  boundary 
layer  structure  can  only  be  accomplished  at  the  expense  of  accurate  predictions  of  the  axial 
Reynolds  stress,  which  can  be  important  in  shock  /  boundary  layer  interactions.  It  is  clear  that 
the  shift  from  LES  to  RANS  brings  about  a  redistribution  of  Reynolds  stresses  that  may  not  be 
correct.  If  the  transition  from  RANS  to  LES  is  designed  so  that  it  occurs  in  the  lower  part  of  the 
logarithmic  region,  then  the  turbulence  intensities  may  be  well-predicted  but  the  overall  mean 
flow  is  accelerated  more  than  it  should  be.  This  leads  to  the  so-called  Tog  law  mismatch’,  which 
plagues  nearly  every  LES/RANS  model. 

A  path  forward  to  resolving  this  issue  is  to  have  a  more  ‘exact’  solution  available  for 
comparison.  To  this  end,  we  have  conducted  wall-resolved  large-eddy  simulations  and  direct 
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numerical  simulations  of  the  Elena  and  Lacharme  [7]  flat-plate  boundary  layer  experiment. 
Resolution  requirements  for  a  wall-resolved  LES  are  much  more  severe  than  for  a  RANS  or 
LES/RANS,  and  the  predictions  are  much  more  sensitive  to  the  level  of  numerical  dissipation 
employed  and  to  the  form  of  the  selected  subgrid  model.  It  also  was  necessary  to  recycle  both 
mean  and  fluctuating  data,  simply  as  there  was  no  RANS  base  state  upon  which  to  super- impose 
recycled  fluctuations.  Several  cases  were  run,  as  summarized  in  Table  9.1.  The  naming 
convention  follows  that  used  in  the  figures  discussed  next.  In  some  cases,  better  resolution  of 
the  flow  in  the  Z  direction  was  expedited  by  reducing  the  Z  extent  of  the  domain  while 
maintaining  a  constant  cell  count.  This  turned  out  to  be  a  bad  idea,  as  it  also  reduced  the 
thickness  of  the  boundary  layer  and  made  some  results  difficult  to  interpret.  The  baseline 
numerical  scheme  used  was  the  LD-PPM  (low  dissipation  piecewise  parabolic  method) 
developed  in  Year  1  of  the  grant.  This  scheme  uses  the  Ducros,  et  al.  vorticity  /  divergence 
switching  function  to  shift  the  inviscid  flux  discretization  from  PPM  in  shock-dominated  or  ffee- 
stream  regions  (switch  =  1)  to  a  fourth-order  central  scheme  in  vorticity-dominated  regions 
(switch  =  0).  The  specific  switch  includes  a  thresholding  constant  that  can  maintain  some  of  the 
PPM  contribution  even  when  the  Ducros  switch  is  nominally  zero.  While  we  have  used  this 
routinely  in  LES/RANS  applications,  it  is  necessary  to  set  this  constant  to  zero  for  LES/DNS 
applications,  as  otherwise,  too  much  numerical  dissipation  can  corrupt  the  solution.  The 
nomenclature  ‘-NT’  refers  to  setting  this  thresholding  constant  (called  ‘cutoff  in  the  Table)  to 
zero. 


Table  9,1:  Cases  run  for  LES/DNS  studies 


Case 

cells  in  X/8 

cells  in  Z/S 

SGS 

cutoff 

discretization 

DNS 

30 

30 

none 

0.1 

LD-PPM 

DNS-NT 

30 

30 

none 

0.0 

LD-PPM 

DNS-PPM 

30 

30 

none 

N/A 

PPM 

DNS- 1/2  dz 

30 

60 

none 

0.0 

LD-PPM 

DNS  fine 

40 

60 

none 

0.0 

LD-PPM 

LES-NT 

30 

30 

Lenormand 

0.0 

LD-PPM 

LES-NT-smag 

30 

30 

Leveque 

0.0 

LD-PPM 

LES-NT-2/3  dz 

30 

45 

Lenormand 

0.0 

LD-PPM 

LES-NT- 1/2  dz 

30 

60 

Lenormand 

0.0 

LD-PPM 

RANS 

30 

30 

N/A 

TVD 

Nomenclature:  DNS  =  direct  numerical  simulation;  LES  =  large-eddy  simulation,  LD-PPM  = 
low-dissipation  piecewise  parabolic  method,  TVD  =  total  variation  diminishing,  cutoff  = 
thresholding  parameter  in  LD-PPM,  dz  =  grid  spacing  in  Z  direction,  SGS  =  subgrid-scale 
closure,  8  =  boundary  layer  thickness;  X  =  streamwise  direction,  Z  =  spanwise  direction 
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Figure  9.1.1  shows  boundary-layer  momentum  thickness  (left)  and  skin  friction  (right).  The 
values  of  the  momentum  thickness  are  smaller  for  the  cases  that  decreased  the  effective  mesh 
spacing  in  the  Z  direction  by  reducing  the  spanwise  extent  of  the  domain.  In  effect,  the  solution 
‘converged’  to  a  thinner  boundary  layer  than  anticipated.  Because  both  the  mean  flow  and  the 
fluctuation  fields  are  rescaled  and  recycled,  there  is  no  direct  way  to  constrain  the  inflow 
boundary  layer  thickness.  None  of  the  LES  or  DNS  solutions  predict  as  high  of  a  skin  friction 


Figure  9.1.1:  Momentum  thickness  (left)  and  skin  friction  (right)  versus  streamwise  distance 


coefficient  as  does  the  RANS  model,  but  as  the  boundary  layer  becomes  thinner,  the  values 
increase.  The  Lenormand  SGS  model  yields  a  higher  skin  friction  level  than  the  Leveque,  et  al. 
[10]  modified  Smagorinsky  model  for  the  same  mesh.  Reducing  the  cutoff  parameter  from  0.1 
to  zero  also  promotes  an  increase  in  skin  friction,  and  the  use  of  the  PPM  method  alone  reduces 
the  skin  friction  greatly,  implying  that  this  scheme  leads  to  excessive  dissipation  of  the  near-wall 
eddy  structures.  The  skin  friction  distributions  also  show  that  there  is  an  ‘adjustment’  period  of 
about  four  boundary  layer  thicknesses  for  the  LES  or  DNS  solutions  before  the  correct  trends 
begin  to  emerge. 

Figure  9.1.2  plots  the  velocity  profile  in  wall  coordinates  (left)  (note  that  the  axes  are  labeled 
incorrectly)  and  in  a  modified  way  in  which  the  velocity  is  normalized  by  the  free-stream 
velocity  and  the  wall  distance  is  normalized  by  the  boundary  layer  thickness  (right)  but  a 
logarithmic  scale  is  still  used.  The  second  way  reduces  the  influence  of  the  general  under¬ 
prediction  in  skin  friction  noted  in  Figure  9.1.1.  Also  shown  in  the  right  figure  as  dashed  lines  is 
the  eddy  viscosity  (dimensional  units)  for  each  of  the  LES  cases. 
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Figure  9.1.2:  Velocity  profile  in  wall  coordinates  (left)  and  in  outer-layer  coordinates  (right) 

There  are  several  points  to  note.  First,  in  the  left  figure,  the  predictions  improve,  relative  to  the 
theoretical  wall  /  wake  law,  as  the  skin  friction  prediction  improves.  The  LES  and  DNS 
predictions  consistently  rise  over  the  wall  /  wake  law  in  the  buffer  region,  while  the  RANS 
prediction  lies  below  the  theory.  The  right  figure  also  includes  Coles’  law  of  the  wall  /  wake, 
transformed  into  outer-layer  coordinates  as  yellow  dots.  Two  such  plots  are  shown  -  one  which 
includes  a  hyperbolic-tangent  blending  of  the  viscous  sub-layer  solution  with  the  logarithmic 
solution  and  one  that  does  not.  The  use  of  outer-layer  scaling  collapses  the  predictions  better. 
This  figure  clearly  shows  that  the  DNS  calculations  with  LD-PPM  and  a  zero  cutoff  and  the  LES 
calculations  with  LD-PPM,  a  zero  cutoff,  and  the  Lenormand  subgrid-scale  model  provide  the 
best  results  in  the  logarithmic  and  wake  regions.  The  Leveque,  et  al.  modified  Smagorinsky 
model  fares  poorly  but  does  lead  to  nearly  vanishing  eddy  viscosity  in  the  laminar  sub-layer.  In 
contrast,  the  Lenormand,  et  al.  SGS  viscosity  does  not  decay  as  rapidly  in  the  buffer  layer  and  as 
a  possible  consequence,  this  model  under-predicts  the  velocity  in  the  laminar  sub-layer. 
Calculations  approach  theory  for  this  model  as  the  mesh  spacing  in  the  Z  direction  is  reduced. 
The  RANS  and  LES/DNS  predictions  in  the  buffer  layer  are  different  and  it  is  not  clear  which 
one  is  correct. 

Predictions  of  Reynolds  shear  stress  and  mean  velocity  (left)  and  rms  axial  and  normal 
fluctuation  velocities  (right)  are  shown  in  Figure  9.1.3  in  outer-layer  coordinates.  Note  that  the 
X-axis  of  the  right  figure  should  also  be  Y 78.  The  actual  boundary  layer  thickness  as  predicted 
by  each  model  is  used  in  the  normalization;  this  reduces  the  scatter  in  the  predictions  for  the 
Reynolds  shear  stress  and  velocity  but  not  for  the  rms  velocities. 
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Figure  9.1.3:  Axial  velocity  and  Reynolds  shear  stress  (right)  and  rms  axial  and  normal  fluctuation 
velocities  (left)  versus  Y/8 


Some  of  the  LES  and  DNS  solutions  agree  better  with  the  theoretical  velocity  profile  than  does 
the  RANS  solution.  The  trends  exhibited  by  the  Reynolds  stress  and  fluctuating  velocity 
components  are  as  expected,  though  there  is  significant  scatter  among  the  predictions. 

The  primary  purpose  of  conducting  these  analyses  is  to  identify  a  ‘reference’  solution  that  can 
then  be  interrogated  to  determine  possible  weaknesses  in  the  LES/RANS  formulation  and 
thereby  focus  potential  improvements.  It  is  unclear  if  any  of  the  obtained  solutions  can  yet  serve 
as  this  ‘truth  model’.  The  fact  that  decreasing  the  Z  extent  of  the  domain  leads  to  a  reduction  in 
the  boundary  layer  thickness  complicates  this  assessment,  as  does  the  fact  that  the  ‘proper’ 
solution  of  the  velocity  in  the  buffer  region  is  not  known.  One  way  of  evaluating  the 
LES/RANS  model  response  may  be  to  compare  forces  exerted  by  the  resolved,  modeled,  and 
molecular  Reynolds  stresses: 

(^)=J((rr)+(rr>+(cl>)'^  (9ii) 

A 

as  they  vary  through  the  boundary  layer.  The  difference  between  ensemble-averaged  forces  from 
a  LES  and  those  from  the  LES/RANS  is  the  form  of  the  modeled  component,  which  for  the 
LES/RANS  methods  considered  herein,  includes  an  unsteady  RANS  contribution,  a  subgrid 
component,  and  a  spatially-varying  blending  function  that  connects  the  two.  The  analogous  LES 
method  includes  only  the  subgrid  component,  so  for  the  forces  exerted  on  the  fluid  for  the  two 
methods  to  be  equal,  the  sum  of  the  resolved  and  modeled  components  must  balance.  In  general, 
this  will  not  occur,  and  the  degree  of  imbalance  represents  model  form  error  as  it  is  reflected  in 
the  time-averaged  solution. 

Based  on  these  statistics  (and  similar  ones  associated  with  heat  and  mass  transfer),  it  should  be 
possible  to  localize  potential  sources  for  model  form  error.  Given  a  particular  LES/RANS 
model  variant,  it  is  possible  to  determine  the  ‘optimal’  model  response  that  leads  to  mean- flow 
equivalence  between  the  LES  and  LES/RANS  solution  on  the  same  mesh.  This  response  may  be 
obtained  by  including  the  difference  between  ensemble-averaged  forces  as  a  source  term  (again 
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with  analogous  expressions  for  the  energy  and  species  equations)  that  steers  the  ensemble 
averages  toward  those  of  the  ‘more  trusted’  LES. 

j ^dtl  +  J (piir «,  +  S„p  - r-  -T™' )-ndA  =  (i- “s)  - (i=™s)  (9. j ,2) 

n  O*  a 

Here,  ^F(LES  ^  will  be  determined  beforehand  from  the  wall-resolved  LES  and  ^LES/RANS  ^  wjH  be 
determined  as  the  calculation  proceeds  by  exponentially-weighted  ensemble-averaging  methods. 
The  final  solution  for^fvLES/RANS^ ,  separated  into  its  components,  represents  the  ‘optimal’  force 

balance  and  thereby  represents  a  focal  point  for  directing  model  improvements.  The  optimal 
LES/RANS  model  will  be  the  one  that  replicates  the  mean-flow  statistics  of  the  wall-resolved 
LES  without  requiring  a  forcing  function. 

10.  Development  of  ‘Data-Mining’  Strategies  for  Analyzing  Reynolds- 
averaged  Navier-Stokes  Closure  Model  Assumptions  Based  on  LES/RANS 
data 

This  last  part  of  the  study  investigated  the  possibility  of  ‘mining’  LES/RANS  solutions  to 
provide  data  that  might  be  used  to  evaluate  and  potentially  improve  RANS-level  models.  As  all 
two-equation  and  second-moment  RANS  models  require  the  solution  of  a  turbulence  dissipation 
rate  equation  (in  some  form),  the  initial  part  of  this  study  details  the  calculation  of  the  turbulence 
dissipation  rate  e  and  related  turbulent  flow  terms  by  using  resolved  flow  data  from  a  numerical 
simulation  utilizing  a  hybrid  LES/RANS  turbulence  closure.  The  estimate  of  the  dissipation  rate 
was  obtained  by  solving  the  Favre-averaged  turbulent  kinetic  energy  equation  throughout  the 
simulated  flow.  The  constituent  terms  of  the  equation,  which  included  fluctuating  and  Favre- 
averaged  flow  variables,  were  obtained  by  ensemble-averaging  various  resolved-scale  flow  data 
over  the  course  of  the  simulation.  Once  the  dissipation  rate  was  obtained,  it  was  used  to  estimate 
the  specific  dissipation  rate,  co,  and  kinematic  eddy  viscosity,  nt,  throughout  the  flow.  These 
were  compared  to  their  fully-modeled  counterparts  obtained  using  Menter’s  BSL  two-equation 
model  in  order  to  analyze  the  accuracy  of  the  latter  in  simulating  complex  flows. 

The  flow  chosen  for  this  analysis  was  a  Mach  2.79  air  flow  over  a  20  degree  wedge  (Figure 
10.1.1.  The  region  at  the  start  of  the  wedge  features  shock-boundary  layer  interaction,  where  the 
accuracy  of  standard  two-equation  turbulence  models  is  questionable.  To  obtain  the  necessary 
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time-averaged  flow  data,  the  simulation  compiled  averages  over  13500  iterations,  representing 
approximately  five  flow-through  times. 


The  Favre-averaged  TKE  equation  reads  as  follows: 


dpk  dpkuj 
dt  dx. 
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Time-averaged  quantities  are  denoted  with  an  overbar,  Favre-averaged  quantities  with  a  tilde, 
and  fluctuating  quantities  with  a  double  prime  symbol,  zv  denotes  Reynolds  stresses  while  tfj 

represents  the  laminar  shear  stress  tensor.  The  individual  fluctuation  and  Favre-averaged 
quantities  were  reconstructed  using  time-averaged  data  as  follows: 
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All  quantities  were  span-averaged  during  post-processing.  Gradients  were  then  computed  using 
a  finite-volume  Green’s  theorem  formulation  at  each  mesh  cell. 


For  this  analysis,  the  molecular  diffusion  term  was  ignored  due  to  the  uncertainty  in  estimating 
the  shear  stress  tensor  very  close  to  the  wall.  The  first  term  on  the  left-hand  side  was  also 
omitted  as  the  analysis  considers  a  quasi-steady  state  solution.  Finally,  because  the  hybrid 
LES/RANS  turbulence  closure  used  in  the  simulation  does  not  resolve  turbulent  eddies  near  the 
wall,  it  was  necessary  to  add  a  modeled  component  to  the  estimate  of  the  Reynolds  stress  tensor 
Ty  when  calculating  the  turbulent  kinetic  energy  production  rate.  The  modeled  stresses  were 

formed  based  on  the  Boussinesq  hypothesis  such  that  the  final  tensor  equations  were: 
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where  //,  and  km  denote  the  eddy  viscosity  and  turbulent  kinetic  energy,  respectively,  obtained 
from  Menter’s  BSL  model. 

The  convection,  production,  diffusion,  pressure  work,  and  pressure  dilatation  were  then 
evaluated  as  described  above  and  used  to  solve  for  the  turbulent  dissipation  rate.  Figure  10.1.2 
shows  the  resulting  values  of  each  calculated  term  in  the  TKE  equation  taken  inside  the 
boundary  layer  at  x  =  -0.2  meters.  The  dominant  term  is  the  turbulent  energy  production  rate, 
which  in  turn  determines  the  dissipation  rate.  The  remaining  terms  tend  to  be  relatively  small, 
with  the  exception  of  the  region  near  the  top  of  the  boundary  layer,  where  the  convection  and 
turbulent  transport  terms  are  briefly  dominant.  However,  these  two  terms  nearly  cancel  out  and 
do  not  influence  the  calculated  dissipation  rate. 


convection 


Figure  10.1.2:  Turbulence  energy  equation  balances  (Mach  2.79  flat-plate 
boundary  layer) 

The  specific  dissipation  rate  and  eddy  viscosity  were  then  estimated  based  on  the  turbulence 
dissipation  rate.  The  specific  dissipation  rate  was  calculated  as  <»  =  £7(0.09  x(£  +  A:m)).  The 
modeled  turbulence  kinetic  energy  km  was  added  to  the  resolved  TKE  value  to  account  for  the 
lack  of  resolved  eddies  near  the  wall,  similarly  to  the  adjustment  made  to  the  resolved  Reynolds 
stress  tensor.  Finally,  the  kinematic  eddy  viscosity  was  given  by  vt  =  (k  +  km )  /  co 

Figures  10.1.3  through  10.1.6  show  progressions  of  the  k,  e ,  co  and  vt  profiles  through  the 
flow,  taken  along  lines  extending  perpendicular  to  the  solid  surface.  Note  that  for  the  TKE  and 
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Figure  10. 1 .3:  Turbulence  kinetic  energy  profiles  throughout  interaction  region 


eddy  viscosity  plots,  the  scale  changes  after  the  contour  at  x  =  -0.04m,  due  to  the  increased  levels 
of  turbulence  caused  by  the  shock  and  recirculation  zone.  The  resolved  values  based  on  the 
above  analysis  are  compared  with  fully  modeled  versions  based  on  a  RANS  simulation  of  the 
same  flow  utilizing  Menter’s  BSL  turbulence  closure. 
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Figure  10.1.4:  Turbulent  dissipation  rate  profiles  throughout  interaction 
region 
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The  resolved  turbulent  kinetic  energy  (Figure  10.1.3)  upstream  of  the  comer  closely  matches  the 
modeled  values  throughout  most  of  the  boundary  layer.  A  significant  deviation  occurs  close  to 
the  wall,  where  the  resolved  TKE  grows  much  larger  as  a  consequence  of  being  augmented  by 


Specific  dissipation 
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Figure  10.1.5:  Specific  turbulent  dissipation  rate  profiles  throughout  the 
interaction  region 

the  modeled  near-wall  values.  Past  the  comer,  the  resolved  TKE  significantly  exceeds  the  model 
within  the  recirculation  zone,  then  approaches  the  modeled  values  again  as  the  flow  re-attaches. 
There  is  also  a  local  spike  in  the  resolved  TKE  values  around  the  shock. 

The  resolved  turbulent  dissipation  rate  (Figure  10.1.4)  is  generally  in  good  agreement  with  the 
modeled  values  in  the  upstream  part  of  the  flow,  although  it  is  difficult  to  make  a  qualitative 
comparison  very  near  the  wall  since  both  values  grow  asymptotically.  Consistent  with  the 
turbulent  kinetic  energy  results,  the  resolved  dissipation  rate  is  relatively  high  inside  the 
recirculation  zone  near  the  comer,  although  the  near-wall  values  at  the  comer  itself  are  lower 
than  the  model  predicts.  The  spike  in  the  resolved  dissipation  rate  around  the  shock  is  more 
pronounced  due  to  the  large  contribution  of  the  pressure  work  in  this  region. 

The  resolved  specific  turbulence  dissipation  rate  (Figure  10.1.5)  follows  a  very  similar  pattern; 
the  main  difference  in  the  upstream  region  of  the  flow  appears  to  be  a  relatively  low  near-wall 
value  of  the  resolved  specific  dissipation,  due  to  a  probable  over-estimate  of  the  TKE.  This 
effect  is  also  evident  past  the  comer,  where  the  modeled  specific  dissipation  rate  is  higher  near 
the  wall  than  the  resolved.  The  erratic  fluctuation  in  the  resolved  data  beyond  the  top  of  the 
boundary  layer  in  the  first  station  is  a  result  of  dividing  by  very  small  TKE  values  and  is  not 
physically  significant. 

The  resolved  kinematic  eddy  viscosity  (Figure  10.1.6)  agrees  with  the  peak  modeled  value  in  the 
upstream  region,  but  is  more  constant  throughout  the  entire  boundary  layer.  The  higher  than 
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Figure  10.1.6:  Kinematic  eddy  viscosity  profiles  throughout  the  interaction  region 

modeled  values  near  the  wall  can  be  explained  by  the  overly  high  resolved  TKE  values; 
however,  it  is  unclear  why  the  resolved  eddy  viscosity  near  the  top  edge  of  the  boundary  layer  is 
significantly  higher  than  the  modeled  form.  Past  the  comer,  the  resolved  eddy  viscosity  tracks 
the  same  pattern  shown  by  the  resolved  TKE,  with  very  high  values  in  the  recirculation  zone. 
The  peak  eddy  viscosity  values  occur  immediately  downstream  of  the  shock,  although  unlike  the 
TKE  and  dissipation  rates,  the  values  upstream  of  the  shock  remain  consistently  minimal. 

These  results  illustrate  the  potential  for  utilizing  LES/RANS  data  to  extract  turbulence 
information  that  could  be  used  in  RANS-level  models.  The  results  are  sensitive  to  the 
interpretation  of  the  RANS  component  within  the  LES/RANS  framework  and  to  the  blending 
between  the  two  methodologies.  Future  work  should  focus  on  ensuring  consistency  between 
these  aspects  of  the  model. 

11.  Summary 

This  research  has  developed  a  new  hybrid  large-eddy  /Reynolds-averaged  Navier-Stokes 
turbulence  closure  strategy  specifically  designed  for  strongly  interacting,  wall-bounded  flows. 
The  model  differs  from  its  predecessor  in  that  the  need  to  pre-calibrate  a  model  constant  is 
removed  through  the  use  of  ensemble-averaged  turbulence  information  to  estimate  an  outer-layer 
turbulence  length  scale.  The  model  has  been  applied  to  a  variety  of  shock  /  boundary  layer 
interactions  and  has  shown  a  good  level  of  predictive  capability  for  both  mean  and  second- 
moment  quantities.  A  specific  result  of  the  shock  /  boundary  layer  interaction  study  is  a  strong 
correlation  between  the  most  probable  time  of  a  fluid  within  the  recirculation  region  formed 
through  shock  interaction  and  the  dominant  low-frequency  signal  of  the  interaction.  This 
provides  evidence  that  the  appearance  of  a  low-frequency  mode  of  separation-shock  unsteadiness 
is  intimately  connected  with  the  structure  of  the  backflow  region  and  the  mean  entrainment 
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patterns.  With  this  knowledge  in  place,  it  may  be  possible  to  predict  low-frequency  dynamics  of 
complicated  interactions  by  examination  of  the  mean  structure  of  the  interactions.  The 
LES/RANS  model  was  also  tested  for  turbulent  flow  over  an  airfoil  near  static  stall  as  an  initial 
step  toward  its  use  in  predicting  dynamic  stall. 

12.  Other  Information 

12.1.  Students  supported  and  degrees  received 

1.  Daniel  Gieseking  (  May,  2009  -  December,  2011)  M.S.  Thesis  “A  New  Hybrid  Large-Eddy 
Simulation  /  Reynolds  Averaged  Navier-Stokes  Model  for  Compressible  Flows”,  currently 
employed  by  Honda  Aircraft. 

2.  Ilya  A.  Zilberter  (January,  2012  -  December,  2012)  Worked  on  data-mining  strategies  using 
LES/RANS  solutions.  Currently  Ph.D  candidate  at  NCSU 

3.  Jianghua  Ke  (December,  2013  -August,  2013)  Worked  on  airfoil  simulations  using 
LES/RANS  for  static  and  dynamic  stall  problems.  Currently  Ph.D  candidate  at  NCSU  supported 
by  new  ARO  grant  with  Dr.  Gopalarathnam  as  PI. 

12.2.  Publications  directly  related  to  this  work. 

Ke,  J.,  Edwards,  J.R.,  “RANS  and  LES/RANS  Simulation  of  Airfoils  under  Static  and  Dynamic 
Stall”  AIAA  Paper  2013-0955,  January,  2013 

Gieseking,  D.  and  Edwards,  J.R.  “Simulations  of  a  Mach  3  Compression-Ramp  Interaction  using 
LES/RANS  Models”  AIAA  Journal,  Vol.  50,  No.  10,  2012,  pp.  2057-2068. 

Gieseking,  D.A.  “A  New  Hybrid  Large-Eddy  Simulation  /  Reynolds-averaged  Navier-Stokes 
Model  for  Compressible  Flows”,  M.S.  Thesis,  Aerospace  Engineering,  NCSU,  December,  2011. 

Gieseking,  D.A.,  Choi,  J.-I.,  Edwards,  J.R.,  and  Hassan,  H.A.  “Compressible  Flow  Predictions 
Using  Improved  LES/RANS  Models”  AIAA  Journal,  Vol.  49,  No.  10,  201 1,  pp.  2194-2209 

Gieseking,  D.A.,  Edwards,  J.R,  and  Choi,  J.-I.  “Simulation  of  a  Mach  3  24-Degree 
Compression-Ramp  Interaction  using  LES/RANS  Models”  AIAA  Paper  2011-5541,  August, 
2011. 

Gieseking,  D.A.,  Choi,  J.-I.,  Edwards,  J.R.,  and  Hassan,  H.A.  “Simulation  of  Shock  / 

Boundary  Layer  Interactions  Using  Improved  LES/RANS  Models”  AIAA  Paper  2010-110, 
January,  2010. 

12.3.  Technology  Transfer 

1.  LES/RANS  methodology  has  been  implemented  into  NASA’s  VULCAN  code  by  Dr.  Robert 
Baurle. 
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2.  LES/RANS  methodology  has  been  used  by  AFRL  to  simulate  reactant  mixing  as  influenced 
by  impinging  shock  waves.  This  work  was  performed  by  Dr.  John  Boles  of  Taitech, 
Incorporated. 

3.  LES/RANS  methodology  has  been  to  compute  reactive  flows  in  the  University  of  Virginia’s 
Scramjet  Combustion  Facility  as  part  of  Dr.  Edwards’  role  as  CFD  lead  in  the  National 
Center  for  Hypersonic  Combined  Cycle  Propulsion. 

4.  LES/RANS  methodology  has  been  implemented  into  a  version  of  Aerosoft’s  GASP  flow 
solver  by  Reece  Neal. 

12.4.  Other  Connections 

1.  Daniel  Gieseking  defended  his  M.S.  thesis  on  December  4,  2011.  He  is  employed  by  Honda 
Aircraft. 

2.  Dr.  Edwards  and  Dr.  A.  Gopalarathnam  submitted  a  successful  proposal  to  ARO  that  will 
develop  new  theoretical  approaches  for  predicting  forces  induced  during  dynamic  stall.  Data 
generated  from  RANS  and  LES/RANS  simulations  will  be  used  to  inform  the  low-order  theory. 

3.  Dr.  Edwards  participated  in  the  1st,  2nd,  3rd,  4th,  and  5th  Shock  Wave  /  Boundary  Layer 
Interaction  Workshop  (2009-2012),  jointly  organized  by  NASA  Glenn  Research  Center  and 
AFRL.  Each  year,  some  aspect  of  the  ARO-sponsored  work  was  presented. 

4.  The  NCSU  team  participated  in  a  blind  study  sponsored  by  AFRL  in  January,  2010.  We 
computed  a  3-D  shock  /  boundary  layer  interaction  experimentally  mapped  at  the  University  of 
Michigan  using  the  ARO-supported  LES/RANS  model  enhanced  by  a  multi-wall  recycling  / 
rescaling  procedure  developed  using  NASA  support.  Our  predictions  were  in  reasonably  good 
accord  with  experimental  measurements. 
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